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1.  SUMMARY 


This  document  reports  on  a  two-year  research  project  by  Kestrel  Technology  and  Honeywell 
Aerospace  Advanced  Technology  to  combine  model-based  development  of  complex  avionics 
control  software  with  static  analysis  of  the  generated  code  to  achieve  assurance  levels  not 
available  from  either  technique  practiced  separately.  We  concentrated  on  two  classes  of 
numerical  algorithms,  linear  digital  filters  (LDFs)  and  integrating  accumulators,  modifying 
existing  versions  of  Honeywell’s  HiLiTE  model-based  development  system  and  Kestrel’s 
CodeHawk  abstract  interpretation  system  to  share  domain  specific  information  about 
implementations  of  these  algorithms.  This  allowed  CodeHawk  to  exploit  model-level 
specifications  and  theoretical  input  bounds  from  HiLiTE  concerning  the  generated  C  code, 
producing  a  much  more  precise  over-approximation  of  the  output  bounds  and  accumulated 
floating-point  error  bounds  than  would  be  possible  with  generic  abstract  interpretation 
techniques.  These  static  analysis  results  were  then  fed  back  into  HiLiTE  to  be  further  exploited 
in  the  formal  verification  of  the  generated  code. 

We  made  some  unexpected  discoveries  during  the  first  year  analysis  of  LDFs  concerning 
analytic  solutions  to  output  and  error  bounds  for  filters  that  changed  the  architecture  of  our  first 
year  approach,  and  redefined  our  focus  for  the  second  year  effort.  In  particular,  Anca  Browne’s 
discovery  of  an  analytic  technique  for  computing  filter  bounds  obviated  the  need  for  the  iterative 
symbolic  evaluation  of  classical  abstract  interpretation,  making  the  time  to  analyze  a  filter 
negligible.  This  contrasts  with  comparable  filter  analyses  performed  for  Airbus  that  take  on  the 
order  of  20  hours  to  get  less  precise  results.  Extending  this  technique  to  floating  point  error 
analysis,  she  discovered  how  to  generate  a  closed-form  solution  for  the  error  as  a  function  of  the 
number  of  iterations  of  the  recurrence  that  are  performed.  This  resulted  in  CodeHawk’s  “error 
bounds”  for  filters  being  reported  back  to  HiLiTE  not  as  an  instance-specific  numeric  interval,  as 
was  originally  envisioned,  but  as  parameters  to  a  generalized  error  function,  allowing  HiLiTE  to 
analyze  multiple  scenarios  with  different  durations  from  a  single  analysis  result.  When  these 
techniques  were  generalized  to  integrating  accumulators  in  year  2,  CodeHawk  was  able  to  report 
both  output  and  error  ranges  as  parameters  to  a  generalized  function. 

We  built  specialized  versions  of  both  HiLiTE  and  CodeHawk  to  exploit  the  shared  model-level 
information  and  exchange  their  results  over  a  JSON  file  interchange  architecture.  This  combined 
system  was  used  first  to  experiment  with  the  joint  modeling  and  analysis  techniques,  and  then  to 
run  a  series  of  test  suites  for  filters  and  integrating  accumulators  that  explored  the  range  of 
variety,  and  analysis  results,  among  those  algorithms.  These  test  suite  results  are  summarized  in 
detail  in  Sections  4.2  and  4.3  below. 
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2.  INTRODUCTION 


The  objective  of  this  project  was  to  develop  a  software  engineering  approach  and  tool  set  that 
combines  the  relative  strengths  of  model-based  and  code-based  analyses  to  provide  a 
comprehensive  and  scalable  analysis  capability  for  numerical  algorithms  in  complex  systems 
software.  Our  approach  was  to  utilize  static  analysis  of  numerical  algorithms  in  combination  with 
model-based  design  technology  to  yield  considerable  benefits  compared  to  applying  static 
analysis  in  isolation.  The  results  of  this  effort  are  intended  to  help  detect  and  mitigate  a  large 
class  of  defects  that  occur  due  to  differences  between  the  intended  semantics  of  design  models 
and  the  actual  behavior  of  the  software. 

Large  and  complex  systems  control  software  increasingly  is  being  developed  using  model-based 
development  tools.  In  this  approach  a  design  model  of  control  algorithms  is  constructed  using 
control  function  blocks  in  a  data  flow  notation.  The  transfer  function  semantics  is  available  at  the 
model  level  along  with  the  semantics  of  feedback  loops,  periodic  rates,  and  design  patterns  such 
as  counters,  validity  status  of  signals,  reset  semantics,  mid-value  selection,  etc.  Some  model- 
based  toolsets  perform  extensive  analysis  of  models  for  early  detection  of  several  types  of  design 
defects  and  have  been  used  as  DO-178B  qualified  tools  in  certification  of  large  scale  commercial 
avionics  systems.  These  tools  can  reason  with  full  semantics  of  control  transfer  functions  and 
discrete  constructs,  but  the  bounds  on  the  ranges  and  errors  can  be  conservative  and  cannot  take 
into  account  the  actual  source  code  constructs  and  numeric  operations  and  the  artifacts 
introduced  by  the  translation  of  the  design  model  into  source  code  and  the  impact  of  object  code 
execution  on  the  processing  architecture.  This  challenge  is  commonly  faced  when  abstractions, 
assumptions,  and  restrictions  are  utilized  to  estimate  the  behavior  of  object  code  at  the  design 
model  level. 

Abstract  interpretation  is  a  fully  automated  mathematical  process  for  discovering  properties  of 
the  execution  behavior  of  a  program.  Because  they  derive  from  a  formal  representation  of  a 
program’s  behaviors,  the  results  of  this  kind  of  static  analysis  have  the  status  of  mathematical 
proofs.  Informally,  abstract  interpretation  operates  by  computing  an  envelope  of  all  possible 
values  of  the  variables  at  each  point  of  the  program.  For  example,  the  possible  values  of  scalar 
variables  can  be  represented  by  a  collection  of  intervals  (one  for  each  variable)  or  a  convex 
polyhedron  (each  dimension  of  the  affine  space  representing  a  program  variable).  These 
envelopes  are  then  employed  as  lemmas  to  prove  safety  properties  of  the  program,  like  the 
absence  of  arithmetic  overflows  or  out-of-bounds  array  accesses.  Some  abstractions  yield  more 
accurate  envelopes,  which  allow  more  safety  conditions  to  be  discharged,  but  may  require  higher 
computational  times.  For  example,  using  intervals  to  compute  the  range  of  variables  is  extremely 
fast  and  they  can  be  applied  to  analyze  million  lines  of  code.  Using  convex  polyhedra  usually 
yields  much  tighter  bounds,  but  the  exponential  complexity  of  the  underlying  algorithms  makes 
this  approach  impractical  for  more  than  a  few  hundred  lines  of  code.  The  complexity  of 
engineering  a  static  analyzer  that  can  analyze  real  codes  effectively  lies  in  the  combination  of 
different  levels  of  abstraction,  so  that  precise  and  costly  abstractions  (like  convex  polyhedra)  are 
only  applied  to  those  parts  of  the  program  that  require  it,  whereas  rougher  but  faster  abstractions 
are  employed  elsewhere.  There  is  no  automated  way  of  doing  this  combination,  which  is  why 
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industrial-strength  static  analyzers  are  usually  specialized  for  a  particular  application  or  family  of 
applications. 

Our  specific  approach  in  this  project  was  to  specialize  Kestrel  Technology’s  CodeHawk  abstract 
interpretation  technology  to  particular  classes  of  numerical  algorithms  in  avionics  systems 
control  software  generated  by  Honeywell’s  HiLiTE  model-based  development  software.  This 
marriage  of  a  model-level  and  a  code-level  tools  allows  both  the  dual  verification  of  model  and  C 
code  correspondences,  as  well  as  considerable  improvements  to  the  precision  of  CodeHawk’s 
abstract  interpretation  due  to  its  exploitation  of  domain-specific,  model-level  information  from 
HiLiTE.  We  built  a  custom  version  of  both  tools  to  exploit  a  round-trip  JSON  integration  from 
model-analysis  to  code,  to  code-analysis  informed  by  model-analysis,  and  back  to  model- 
analysis  confirmed  by  code-analysis. 

2.1.  Original  SOW  Vision 

The  original  SOW  envisioned  specializing  this  combined  model/code  analysis  architecture  to  a 
specific  class  of  numerical  algorithms,  Linear  Digital  Filters  (LDFs),  in  the  first  year,  then  using 
the  results  of  this  experience  to  explore  other  classes  of  numerical  algorithms  over  a  specific 
abstract  domain,  zonotopes,  in  the  second  year.  Our  first  year  exploration  of  LDFs,  however,  led 
to  some  unexpected  discoveries  about  the  abstract  interpretation  of  filters,  when  the  abstract 
interpreter  knows  (via  model-level  input)  that  the  C  code  it  is  analyzing  actually  implements  a 
filter  function.  This  led  to  a  radical  redesign  of  the  analyzer,  and  a  shift  in  focus  for  the  abstract 
domains  first  envisioned  in  the  SOW.  This  re-orientation  also  led  us  to  re-conceptualize  the 
second  year  effort  toward  a  particular  class  of  numerical  algorithms  (integrating  accumulators), 
instead  of  an  abstract  interpretation  domain  (zonotopes). 

2.2.  Discoveries  about  Abstract  Interpretation  of  Filters 

Abstract  interpretation  attempts  to  compute  an  approximation  of  a  program’s  behaviors  over  all 
possible  execution  paths,  without  any  actual  inputs,  using  only  the  semantics  of  the  programming 
language  to  inform  its  analysis.  Without  knowledge  of  inputs,  or  of  what  algorithms  the  program 
is  executing,  the  inferred  behaviors  will  necessarily  be  a  conservative  over-approximation  of  the 
actual  behaviors.  This  over-approximation  can  often  be  used  to  prove  properties  about  the  actual 
behaviors  of  the  program,  however,  if  the  properties  are  also  true  of  the  over-approximation  of 
those  behaviors.  We  recast  this  problem  for  LDFs  by  using  HiLiTE  to  1)  inform  CodeHawk  that 
the  C  code  it  is  analyzing  implements  a  LDF,  and  2)  provide  theoretical  bounds  on  the  values  of 
all  inputs  to  the  filter.  CodeHawk  is  then  able  to  compute  an  approximation  of  actual  bounds  of 
the  output  values  that  is  much  more  precise  that  would  be  possible  without  these  hints. 

2.2.1.  Filter-specific  widening 

Filters  are  a  special  class  of  iterating  numerical  algorithms  that  use  the  results  of  prior  iteration 
steps  to  constrain  the  computation  of  the  next  step.  This  leads  to  a  natural  tamping  down  of  the 
input  variety  and  causes  the  bounds  of  possible  outputs  to  converge  to  a  smaller  range.  Without 
this  special  knowledge,  uninformed  abstract  interpretation  will  continue  to  “widen”  the  iterated 
values  in  search  of  a  much  more  conservative  convergence.  The  challenge  for  abstract 
interpretation  of  filters  is  in  developing  a  special  widening  analysis  for  filters  that  takes  account 
of  the  unique  cancellation  effects  that  one  output  has  on  the  next  input. 
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We  originally  anticipated  that  our  filter-specific  analyzer  would  be  specialized  to  particular 
classes  of  filters  (exploiting  a  filter-class  parameter  from  HiLiTE)  so  that  widening  could  exploit 
the  unique  cancellation  characteristics  of  these  classes.  We  discovered,  however,  that  this  class- 
specialization  is  unnecessary.  We  developed  an  incremental  widening  algorithm,  taking 
cancellation  into  account,  which  appears  to  work  well  for  all  linear  filters.  Moreover,  we 
discovered  that  a  stable  bound  on  the  filter  output  can  be  computed  symbolically,  via  closed  form 
solution,  thus  obviating  the  need  for  the  traditional  iterative  widening  of  classical  abstract 
interpretation. 

2.2.2.  Floating  Point  Error  Bounds 

Another  original  goal  of  this  project  was  to  take  account  of  floating-point  round-off  errors  when 
real-numbered  models  are  implemented  as  floating-point  code.  Theoretical  bounds  on  filter 
outputs  derived  from  models  over  the  real  numbers  will  often  not  anticipate  possible  round-off 
effects  in  the  generated  C  code,  which  will  vary  with  machine  precision.  Abstract  interpretation 
typically  performs  its  analysis  over  the  infinite  precision  reals  as  well,  so  a  combination  of 
model-level  and  code-level  analyzers  might  still  miss  this  important  delta. 

We  had  originally  planned  to  run  the  analysis  for  each  filter  against  both  CodeHawk’s  native 
real/rational  abstract  domain  and  a  newly  implemented  floating-point  abstract  domain,  and  then 
compare  the  output  range  approximations  to  compute  an  approximation  of  the  round-off  error 
exposure  for  floats.  Our  initial  implementation  of  this  comparison  revealed  that  the  accumulated 
floating-point  error  range  is  a  function  of  the  number  of  sampling  periods  for  the  filter  and  thus 
could  be  computed  independently.  When  this  float-error  accumulation  is  combined  with  the 
cancellation  effects  of  subsequent  inputs,  there  is  an  inherent  stabilization  of  the  error 
accumulation.  This  result  turned  out  to  be  surprising  in  that  the  Honeywell  team  had  estimated 
error  accumulations  to  be  in  the  1  %- 1 0%  range  per  1,000,000  samplings.  The  analysis  results 
yielded  error  ranges  that  were  orders  of  magnitude  lower. 

It  appears  that  filters’  inherent  functionality  in  stabilizing  sensor  output  works  to  stabilize,  rather 
than  accumulate,  round-off  errors.  We  were  further  able  to  derive  an  optimal  function  to  compute 
the  error-accumulation  for  a  filter-instance/sampling-period  combination,  so  rather  than  reporting 
error  bounds,  CodeHawk  now  reports  the  constant  terms  to  this  function,  specific  to  the  filter 
instance  (input  bounds)  under  analysis,  rather  than  an  instance-specific  range.  The  Honeywell 
team  can  now  experiment  with  various  sampling  periods  for  this  filter  instance  off-line  without 
having  to  rerun  the  CodeHawk  analysis. 

2.3.  Re-orientation  of  Year  2 

The  original  idea  in  the  SOW  for  the  second  year  of  this  project  was  to  extend  the  LDF  analysis 
to  an  interestingly  different  class  of  numerical  algorithms.  This  was  originally  characterized,  not 
by  an  algorithm  class  (like  LDFs),  but  by  a  new  analysis  domain:  zono topes.  But  our  first  year 
experience  concluded  that  this  work  would  be  more  valuable  to  the  aerospace  community  if  the 
focus  was  on  some  particular  numerical  algorithm  class  that  has  practical  interest  for  aerospace, 
rather  than  on  a  particular  abstract  interpretation  domain.  Our  somewhat  surprising  discoveries 
about  floating-point  error  convergence  in  FDFs  led  to  an  interest  in  the  converse  class  of 
numerical  algorithms  where  such  errors  do  accumulate,  often  undetected  at  design  time.  We 
started  with  an  “accumulating  integrator”  class  of  algorithms  that  was  the  source  of  a  floating- 
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point  error  design  flaw  in  the  Patriot  Missile  control  software  that  failed  to  detect  an  incoming 
Iraqi  Scud,  allowing  it  hit  the  Army  barracks  in  Dhahran,  Saudi  Arabia  in  1991  (killing  28 
Americans).  If  there  had  been  an  abstract  interpreter,  specialized  to  this  algorithm  class  over  the 
“arithmetic-geometric  progression”  domain,  the  floating-point  round-off  error  could  have  been 
discovered  at  design/implementation  time,  before  the  software  was  deployed. 

We  began  with  a  particular  case  study:  the  Patriot  Missile  failure,  then  generalized  to  other 
numerical  algorithms  in  this  general  class.  The  new  abstract  domain  is  the  arithmetic-geometric 
progression  domain.  It  is  useful  for  analyzing  the  accumulation  of  round-off  errors  in  algorithms 
that,  unlike  filters,  accumulate  rather  than  stabilize  their  errors.  We  got  a  head  start  on  the 
analysis  of  this  class  of  algorithms  through  characteristics  that  they  share  with  the  analytic 
solutions  discovered  for  LDFs.  Although  the  floating-point  errors  accumulate  for  this  class 
(leading  to  application  failures  if  the  bounds  are  not  correctly  anticipated),  rather  than  converge 
as  they  do  for  LDFs,  the  error-bounds  can  still  be  characterized  by  a  single,  analytic  function 
over  the  number  of  iterations  of  the  accumulator.  This  allowed  us  to  report  the  error-bounds 
output  for  accumulators  in  the  same  format  as  for  our  LDF  work.  We  further  discovered  that  the 
range  bounds  on  the  accumulator  output  values  could  be  more  generally  represented  as 
coefficients  to  an  analytic  range  function,  rather  than  a  specific  numeric  interval,  so  the  code¬ 
analysis  results  passed  from  CodeHawk  back  to  HiLiTE  are  more  general  for  accumulators  than 
for  LDFs. 
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3.  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 


The  goal  of  this  project  was  to  build  a  combined  system  of  model-level  analysis  and  code-level 
analysis  that  would  allow  each  analyzer  to  exploit  artifacts  from  the  other.  Since  we  began  with 
two  existing,  unrelated  model  analysis  and  code  analysis  systems,  the  Statement  of  Work  for  the 
project  was  organized  around  three  major  subtasks:  1)  enhancements  to  Kestrel  Technology’s 
CodeHawk  abstract  interpretation  software  to  exploit  model-level  inputs,  2)  enhancements  to 
Honeywell’s  HiLiTE  modeling  and  code  generation  software  to  generate  those  inputs  and  exploit 
the  resulting  analysis,  and  3)  an  integration  architecture  for  allowing  the  enhanced  versions  of 
both  tools  to  work  together  as  a  single  tool,  presenting  a  unified  analysis  to  the  end  user.  The 
following  three  subsections  discuss  the  starting  methods,  assumptions,  and  procedures  for  each 
of  these  subtasks,  respectively. 

3.1.  Control  Algorithm  Capture  in  Models  and  Analysis  /  Code  Verification 

Use  of  Model-Based  Development  (MBD)  techniques  for  software  development  for  control 
algorithms  has  become  a  common  practice  in  avionics  for  systems  of  small  to  large  complexity. 
In  this  approach,  a  design  model  of  control  algorithms  is  constructed  using  control  function 
blocks  in  a  data  flow  notation.  The  control  function  blocks  include  arithmetic  operators,  digital 
filters,  integrators,  timers,  limiters,  etc.  that  are  composed  in  various  patterns  such  as 
proportional-integrator  (PI)  controllers  that  utilize  various  aspects  of  control  including  feedback 
error  compensation,  integrator  anti-windup  logic,  and  shaping  the  characteristics  of  the  inputs 
from  sensors/plant  and  the  command  to  actuators  to  meet  the  control  requirements.  Besides  the 
control-related  constructs,  the  controller  software  also  includes  discrete,  time-dependent  logic 
including  timers  and  counters  for  implementing  mode  and  states  changes. 

The  controller  design  is  captured  in  design  models  using  tools  such  as  MATLAB  Simulink  and 
SCADE.  Source  code  in  C-language  is  generated  from  these  models  that  is  then  used  as  the 
airborne  code  in  avionics  systems.  Processes  and  automation  tools  have  been  deployed  to  certify 
these  systems  to  the  DO-178B  (and  the  revised  DO-178C)  standard.  The  processes  follow  the 
verification  objectives  defined  in  DO-178C  that  start  with  the  software  high-level  requirements 
being  refined  into  design  models  from  which  source  code  is  generated  that  gets  compiled  and 
linked  into  airborne  executable  object  code.  Figure  1  illustrates  this  process.  Note  that  there  are 
several  verification  objectives  defined  in  DO-178C  for  the  development  artifacts  -  only  a  few  of 
the  objectives  are  shown  in  this  figure  for  the  design  models  and  the  source  code. 
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DO-178C:  A-4.2,  3,  5  DO-178C:  A-5.6 


Figure  1.  Software  Development  Process  and  Verification  Objectives 

The  objective  shown  in  Figure  1  for  the  design  model  is  to  verify  that  the  low-level  requirements 
(i.e.,  the  elements  of  the  design  model)  are  accurate,  consistent,  and  verifiable.  This  generally 
includes  absence  of  overflow  (e.g.,  division  by  zero),  untestable  conditions,  and  unreachable 
model  elements.  The  source  code  verification  objectives  are  also  quite  similar,  with  the  addition 
of  memory  safety  specific  to  source  code  constructs,  memory  usage  limits,  etc.  A  point  to  note 
here  that  these  verification  objectives  are  quite  similar  and  are  applied  independently  at  the 
design  model  level  and  source  code  level  -  with  different  methods  and  analysis  tools  at  these  two 
levels  to  perform  the  verification. 

Need  to  achieve  precise  and  consistent  verification  results  across  model-level  and  code-level 
analysis:  At  the  core  of  the  static  analysis  techniques  at  the  model  or  code  level  is  the  derivation 
of  conservative  range  and  numerical  error  bounds  on  the  variables  in  the  algorithm  -  using  these 
bounds  various  properties  can  be  proven.  In  complex  numerical  algorithms,  this  is  a  hard 
problem.  For  example,  for  transfer  functions  such  as  linear  digital  filters,  there  is  a  recurrence 
relation  of  the  output  of  the  filter  to  the  previous  step’s  output  and  input  values,  which  in  turn 
depend  upon  the  values  in  the  step  before  that.  It  turns  out  that  general-purpose  static  analysis 
techniques  give  poor  results  on  linear  digital  filters.  What  is  needed  are  specialized  techniques 
to  derive  sound  and  precise  results.  The  application  of  such  specialized  techniques  needs  to  draw 
upon  the  semantics  of  the  modeling  constructs  so  that  the  tools  can  identify  which  parts  of  the 
design  or  code  a  particular  specialized  analysis  needs  to  be  applied  and  what  are  the  parameters 
of  interest  in  the  particular  instance  of  usage.  This  is  a  common  practice  in  use  of  formal  analysis 
and  model  checking  tools  -  providing  auxiliary  lemmas  and  parts  of  proof  to  the  tool  to  make  the 
problem  tractable. 

A  problem  observed  with  the  current  state  of  the  model-level  and  code-level  analysis  is  that  the 
users  get  different  answers  from  the  tools  applied  at  these  two  levels  due  the  practical  limitations 
of  the  analysis  capabilities  of  the  tools.  For  example,  a  control-theory  based  analysis  may  predict 
reasonable  range  bounds  at  the  output  of  a  filter  but  a  code  analysis  tool  may  derive  very  large 
bounds  for  the  filter  since  that  tool  does  not  apply  a  specialized  analysis  technique  for  this  code. 
This  inconsistency  (and  large  bounds  reported  by  the  code  analysis)  creates  a  problem  for  the 
user  to  understand  and  resolve  the  discrepancy.  Furthermore,  large  bounds  can  lead  to  false 
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alarms  and  trigger  unnecessary  additional  work  steps  in  the  certification  process  to  manually 
analyze  and  resolve  each  situation. 

The  research  goal  of  this  project  is  to  investigate/develop  a  model-based  software  engineering 
approach  and  tool  set  that  combines  the  relative  strengths  of  model-based  and  code-based 
analyses  to  provide  a  comprehensive  and  scalable  analysis  capability  for  numerical  algorithms  in 
complex  systems  software. 

3.1.1.  Use  of  Semantics  at  the  Design  Model  Level  to  Aid  Analysis 

Models  use  control  transfer  functions  such  as  filters,  integrators,  gains  to  build  a  control 
algorithm.  Figure  2  shows  a  simple  Simulink  model  that  is  part  of  a  control  algorithm  feedback 
loop  -  specifically  compensation  of  the  feedback  error  from  the  plant  sensor  to  meet  the  control 
objectives  of  rise  time  and  overshoot.  There  are  two  lead-lag  filters  in  series;  instances  of  the 
generic  lead-lack  filter  block  with  specific  values  of  the  time-constant  parameters. 


Figure  2.  Example  Model  of  Control  Algorithm  using  Lead-Lag  Filters 

The  transfer  function  of  the  lead  lag  filter  is  shown  in  Figure  3.  The  control  algorithm  analysis 
and  simulation  in  closed-loop  mode  starts  with  the  continuous  domain  Laplace  transform.  During 
the  model-based  development  process,  that  is  translated  to  the  transfer  function  in  the  discrete 
domain  (Z-domain)  that  is  implemented  in  the  controller  software.  This  is  a  necessary  step  since 
the  controller  software  is  implemented  as  a  periodic  thread  where  the  transfer  function  is 
evaluated  at  each  periodic  step  (frame).  From  the  Z-domain  transform,  we  derive  a  difference 
equation  that  is  then  implemented  in  the  software.  The  terminology  used  in  Figure  3  is  shown  on 
the  right  hand  side:  u  is  the  input  andy  is  the  output  of  the  filter.  The  subscripts  n,  n-1,  n+1 
denote  the  relative  number  of  a  particular  time  step  in  the  sequence. 
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Continuous  domain: 


u 


Tn  s  +1 


Tds  +1 


7*5+ 1 


Discrete  domain: 


H(-\  _  fc3+M  1 

“  *1+fc2z_1 


leadLag 


yn  -  Current  output 
yn-i  -  Previous  output 
un  -  Current  input 
un_!  -  Previous  input 
Ts  -  Sampling  period 


Coefficient  definitions 

(linearity  assumption:  all  coefficients  are  constantfor  a  filter  instance  in  a  model) 


2Tj 

ki  =  1  +  ~t7 


kj  —  1  —  ■ 


2  Ta 


k3  =  1  +  ■ 


2T„ 


=  1- 


2T„ 


Difference  equation:  yn  =  (fc3un  +  ktun_ A  —  k2yn-f)/k1 


Figure  3.Transfer  Function  for  a  Lead  Lag  Filter  and  Derivation  of  Difference  Equation 

Using  the  Z-domain  transfer  function  H(z),  a  difference  equation  can  be  derived  as  shown  in  the 
figure.  The  difference  equation  computes  the  value  of  y  in  the  current  step  (y„)  in  terms  of  output 
value  in  the  previous  step  and  input  values  in  the  previous  and  current  steps.  The  coefficients  kl 
to  k4  are  defined  in  terms  of  the  sampling  period  Ts  and  the  filter  time  constants:  Tn  is  the 
numerator  time  constant  and  Td  the  denominator  time  constant  of  the  transfer  function.  Note  that 
the  terms  w  and  v,  in  the  difference  equation  correspond  to  variables  in  the  code  where  un-i  and 
Vn-j  are  state  variables  that  hold  the  values  of  u  and  v  from  the  previous  step.  The  difference 
equation  and  the  form  of  its  corresponding  translation  to  code  are  very  important  from  the 
perspective  of  the  code-level  analysis  that  is  described  in  another  section  in  this  document. 

3.1.2.  Exploration  of  Model-Based  Code  Generation  Approach 

Specific  model-level  encapsulation  guidelines  and  code  generation  options  are  required  to  ensure 
the  generated  code  will  not  exhibit  intermingling  of  the  numerical  algorithm’s  code  with  the  code 
from  other  parts  of  the  model.  This  presents  an  unnecessary  difficulty  for  the  code-level  analysis 
since  the  intermingled  code  cannot  be  cleanly  mapped  to  the  appropriate  abstract  domain 
corresponding  to  the  algorithm  characteristics. 


Using  MATLAB  Simulink,  we  perfonned  experiments  with  different  options  of  modeling 
abstractions  and  code  generation  and  selected  a  strategy  that  is  best  suited  for  integrated  static 
analysis  of  models  and  code.  For  each  filter  block,  we  construct  a  reference  model  for  simulation 
of  its  discrete  behavior  and  code  generation.  The  reference  model  is  built  with  basic  arithmetic 
blocks  and  unit  delay  blocks  in  Simulink. 


From  the  continuous  transfer  function  of  a  filter  block,  we  first  derive  its  discrete  transfer 
function  with  parameters  expressed  in  the  continuous  domain.  The  correctness  of  the  discrete 
transfer  function  can  be  verified  by  the  difference  equation  converted  from  it.  There  are  several 
methods  to  transfonn  continuous  transfer  functions  to  discrete  and  vice  versa.  We  use  the 
bilinear  transform  which  essentially  substitutes  “v”  in  the  continuous  transfer  function  using  the 
first  order  approximation: 


2  1  -  z1 
Tt  +  z1 


(1) 
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For  example,  from  the  continuous  transfer  function  H(s)  of  the  lead-lag  filter,  we  have 


H(s)  = 


TnS+1 

Tds+1 


W(z)  = 


k3+k4z 

—  1  f 

ki+k2Z  1 


^2 


l 


2 Ta 
Ts  ’ 


(2) 


Hence  the  difference  equation  is  yn  =  ( fe3u„  +  k4un_1  -  k2yn-i)/k1.  Note  that  T„  and  Td  are  time 
constants  defined  in  the  application  model  for  each  filter  instance  similar  to  the  way  a  constant 
block’s  value  is  specified.  The  sample  time,  Ts,  is  linked  to  the  application  model’s  sample  time 
by  Simulink  for  all  filter  instances. 


Figure  4  shows  a  Simulink  model  that  can  be  referenced  by  application  models.  An  application 
model  employs  a  filter  by  connecting  signals  used  in  an  application  to  a  model  reference  that 
links  to  the  implementation  of  the  filter. 


Referenced  Model 
Program: 

CPU/Process:  C:/TestExample 
DLL:  ModelReferenoe 
Filename:  RefMdl_LeadLagFilter.mdl 

HAM  Library  Version:  7.5.0 
Revision:  - 

Date:  30-0ct-2014  at  14:14:21 
Model  Update  Rate:  Inherited 
[FixedStepDiscrete] 

Trace  Tag:  tbd 
Code  Generation:  Inline 
DLL  Base  Address:  0 
Thread  Name: 

Page:  1  of  1 


Figure  4.  Reference  model  for  lead-lag  filter  block 

The  reference  model  is  constructed  by  the  following  steps: 

1 .  Create  a  new  Simulink  model  file. 

2.  Set  up  the  reference  model  parameters  as  shown  in  Figure  5. 

a.  Navigate  to  “View”-“Model  Explorer”  to  open  the  Model  Explorer  window. 

b.  Click  “Model  Workspace”  in  “Model  Hierarchy”  menu  on  the  left  side  of  the 
window. 

c.  Click  “Add  Simulink  parameter”  button  on  the  top  tool  bar. 

d.  Double  click  the  parameter  to  set  its  name/default  value/type/dimension. . . 
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e.  Add  all  the  parameters  (separated  by  comma)  into  the  model  argument 
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Figure  5.  Setting  reference  model  parameters 

3.  Create  the  Simulink  block  diagram  according  to  the  difference  equation. 

4.  Assign  constant  blocks  with  corresponding  parameters  from  workspace. 


Figure  6  shows  an  application  model  with  lead-lag  filter  blocks.  The  code  for  the  lead-lag  filter  is 
generated  as  a  separate  function.  The  coupling  between  the  main  model’s  code  and  the  filter 
function  is  achieved  via  the  variables  corresponding  to  the  input  and  output  of  the  filter  function. 
Additional  relevant  parameters  are  specified  when  the  model  reference  is  instantiated,  allowing 
the  application  model  to  define  the  time  constants  for  each  instance  and  enforcing  a  consistent 
same  period  between  the  application  model  and  all  filter  instances. 


Lead  Lag  Filter 


Figure  6.  An  Application  Model  Containing  a  Lead-Lag  Filter  Function 

The  code  generated  for  the  model  is  shown  below.  The  input  variables  are  errorcomp  B. Limit  1, 
errorcomp  p .  constant  i  vaiue,  and  Reset.  The  output  variables  is  rtb_cmp3_y.  There  are  two 
State  Variables,  errorcomp_DWork.cmp3_DWORKl.rtb  and  errorcomp_DWork.cmp3_DW0RKl.rtdw  .  The 

sample  time,  errorcomp  p.  cmp3_rtp_rate,  is  automatically  inherited  from  the  application  thread 
level.  Finally  the  filter  instance  time  constants,  errorcomp  p.  cmp3  rtp  Tn  and 
errorcomp  p.  cmp3_rtp  Td,  are  specified  by  the  application  model  designer. 


void  errorcomp  step (void) 

{ 

errorcomp_B . Limit_l  =  saturate_dbl (errorcomp_B . Suml , 

errorcomp  P. Limit  1  LowerSat, 
errorcomp_P . Limit_l_UpperSat )  ; 

/*  call  to  lead  lag  filter  function  for  1st  filter  instance  */ 

RefMdl_ResetLeadLagFilter_p ( Serrorcomp  B . Limitl , 

SerrorcompP. Constant_l_Value,  &Reset, 
&rtb_cmp3_y, 

& (errorcomp_DWork. cmp3_DW0RKl . rtb) , 

& (errorcomp  DWork.cmp3  DW0RK1 . rtdw) , 

errorcomp  P.cmp3_rtp  rate. 
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errorcomp_P . cmp3_rtp_Tn, 
errorcomp_P . cmp3_rtp_Td) ; 

J _ 

It  is  important  that  the  referenced  model  code  closely  reflects  the  difference  equation  for  the 

filter  function.  Other  implementations  may  incorporate  coefficients  into  feedback  loops  to 
produce  the  same  results  with  variable — rather  than  constant — coefficients  that  throw  off  code¬ 
level  analysis. 

Some  filters  (those  named  Variable<filter  type>)  take  tau  (or  Tn  or  Td)  as  an  input  signal.  This 
would  allow  the  model  to  vary  tau  continuously,  which  is  not  supported  by  these  analyses.  The 
coefficients  need  to  be  fixed  for  the  analysis.  In  practice,  however,  tau  is  never  varied 
continuously  since  that  can  make  the  entire  control  analysis  non  linear.  In  general,  the  tau  input 
is  a  way  to  decouple  the  filter  for  more  flexible  models.  For  example,  product  line  modeling 
reuses  models  for  different  variants  of  a  component.  The  architecture  is  the  same  for  all  part 
numbers,  but  details  change  between  parts.  In  such  a  scenario,  a  different  tau  is  specified  for 
different  configurations,  but  it  is  fixed  for  each.  In  other  usages,  where  multiple  modes  exists  in 
the  avionics  subsystems,  tau  can  change  across  mode  changes,  but  is  constant  within  a  single 
mode.  This  allows  separate  piecewise  analyses  to  be  performed  for  each  mode. 

3.2.  Abstract  Interpretation 

Given  a  model  and  code  generated  from  the  model,  we  want  to  provide  mathematical  evidence 
that  the  code  faithfully  implements  the  model.  For  linear  digital  filters  (LDFs),  the  code 
implements  a  linear  recurrence  equation  derived  from  a  continuous-domain  model.  A  filter  takes 
a  stream  of  inputs  (e.g.  a  stream  of  sensor  readings),  and  produces  a  stream  of  output  values  in 
accord  with  the  linear  recurrence.  An  accumulating  integrator  takes  a  stream  of  inputs  and  also 
produces  a  stream  of  outputs,  again  according  to  a  linear  recurrence,  but  the  output  values  and 
their  floating  point  errors  can  grow  without  bound  with  the  number  of  iterations. 

Two  main  analysis  issues  arise:  (1)  what  is  the  range  of  possible  output  values  considered  over 
all  possible  input  streams,  and  (2)  what  is  the  range  of  possible  errors  in  the  output  stream  due  to 
the  floating  point  implementation  of  real-valued  variables  in  the  continuous  model? 

Generally,  static  analysis  reasons  about  all  possible  behaviors  of  a  program.  Since  there  are 
usually  a  very  large  (or  even  infinite)  set  of  possible  inputs,  it  is  usually  not  feasible  to  reason 
about  all  the  inputs  individually.  Instead,  static  analysis  treats  all  inputs  at  once  by  characterizing 
their  common  structure,  via  type  information  and  logical  constraints.  The  goal  is  to  symbolically 
simulate  the  execution  of  the  program  over  this  input  characterization,  yielding  a  characterization 
of  all  possible  behaviors.  Typically,  programs  have  a  very  complex  set  of  behaviors  and  so  it  is 
only  feasible  to  generate  an  approximation,  analogous  to  finding  the  hull  of  a  complex  set  of 
points  in  space.  If  the  approximation  is  an  over-approximation  (i.e.  an  upper  bound),  then  any 
property  that  we  can  prove  about  the  over-approximation  also  holds  for  all  behaviors  of  the 
program. 

The  general  theory  of  static  analysis,  called  Abstract  Interpretation,  was  developed  in  the  1970's 
by  Cousot  and  Cousot  [1].  The  technique  is  widely  used  and  can  be  efficient  enough  for 
commercial  use.  One  of  the  success  stories  of  abstract  interpretation  is  the  complete  analysis  of 
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the  Airbus  380  flight  control  software  for  out-of-bound  memory  access  and  arithmetic  exceptions 
[2],  CodeHawk  is  a  state-of-the-art  implementation  of  abstract  interpretation. 

The  key  idea  of  abstract  interpretation  is  abstract  domains  that  allow  computing  an  over¬ 
approximation  of  the  behaviors  of  a  given  program.  An  abstract  domain  provides  (1)  an  abstract 
type  to  represent  concrete  program  states,  and  (2)  abstract  functions  to  represent  the  effect  of 
concrete  state-changing  actions.  Rather  than  simulate  the  concrete  program,  abstract 
interpretation  uses  abstract  domains  to  construct  and  simulate  an  abstract  version  of  the  program. 
The  simulation  is  computationally  cheaper  in  the  abstract  program  because  the  abstraction 
throws  away  infonnation.  However  the  very  act  of  throwing  away  infonnation  causes  a  loss  of 
precision  in  the  simulation/analysis  process.  Hence  there  is  a  tradeoff  in  abstract  interpretation 
between  computational  complexity  of  the  analysis  and  precision  of  the  results. 

To  illustrate  the  loss  of  precision  due  to  abstraction,  suppose  that  we  want  to  analyze  the 
following  program  fragment,  where  input  x  ranges  over  the  integers  {1,2,3, 4,5}: 

y  =  3*x; 

if  y  >  20  then  ... 

If  we  abstract  each  variable  to  an  interval  over  the  integers,  then  we  can  precisely  represent  x  in 
the  initial  and  final  state  of  this  simple  program  by  the  interval  [1,5]  which  gives  a  minimum  and 
maximum  value  that  the  variable  can  take  on  in  any  execution  of  the  program.  However,  to 
reflect  the  action  of  multiplying  by  3,  the  abstract  domain  would  multiply  the  interval  by  3, 
resulting  in  the  interval  [3,15]  for  y.  Clearly,  we  have  lost  infonnation  in  this  abstraction,  since 
we  can  only  represent  the  possible  values  of  y  via  a  simple  interval  (rather  than  the  precise  set 
{3,6,9,12,15}).  On  the  other  hand,  the  abstraction  does  allow  us  to  cheaply  compute  some  kinds 
of  information  about  the  concrete  program.  In  the  example,  we  can  symbolically  evaluate  the 
condition  to  false  in  the  abstract  domain,  so  we  know  that  the  then-branch  can  never  be  executed. 

An  abstract  interpretation  engine  comes  with  a  library  of  standard  abstract  domains  that  can  be 
used  to  analyze  a  broad  range  of  problems.  The  most  common  abstract  domain  is  numeric 
intervals,  as  illustrated  above,  which  abstract  the  possible  values  of  a  numeric  variable  at  a 
program  control  point  by  a  bounding  interval.  Another  common  abstract  domain  uses  a  set  of 
linear  constraints  (i.e.  an  enclosing  polyhedron)  to  over-approximate  the  joint  values  of  several 
variables.  The  interval  domain  is  computationally  cheap  and  scales  well  to  very  large  programs, 
although  it  can  become  imprecise.  In  contrast,  a  polyhedral  abstract  domain  provides  good 
precision  since  it  relates  several  variables,  but  can  be  very  expensive  to  apply. 

An  abstract  domain  is  used  to  generate  a  specific  kind  of  invariant  at  each  program  control  point 
(e.g.  the  control  point  between  assignments  in  C).  The  interval  domain  infers  interval  invariants 
for  each  program  variable  at  each  control  point.  The  polyhedral  domain  infers  a  polyhedral 
envelope  for  some  of  the  program  variables  at  each  control  point. 

To  analyze  a  given  program,  an  abstract  interpreter  starts  with  an  abstract  representation  of  the 
input  and  forward  simulates  that  representation  through  the  abstracted  actions  of  the  program. 

The  simulation  continues  until  a  fixpoint  is  reached,  which  is  expressed  as  a  sound  invariant  at 
each  control  point  in  the  program.  The  presence  of  loops  complicates  the  analysis,  since  it  may 
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be  necessary  to  simulate  around  a  loop  an  unlimited  number  of  times.  Abstract  interpretation 
uses  a  widening  operator  to  speed  up  convergence  of  the  simulation  to  a  fixpoint,  by  generalizing 
the  current  abstract  representations.  The  result  will  generally  be  a  fixpoint  (set  of  invariants  at 
control  points)  that  is  not  a  least  fixpoint  (the  strongest  expressible  set  of  invariants  at  control 
points  of  the  program). 

Part  of  the  art  of  abstract  interpretation  is  carefully  choosing  where  in  a  program  to  apply  which 
abstract  domain  -  where  precision  is  needed  (as  in  an  inner  loop)  we  might  apply  a  polyhedral 
domain,  but  elsewhere  use  intervals.  The  goal  is  to  gain  as  much  information  as  possible  about 
the  program  while  minimizing  analysis  time  and  space.  Another  part  of  the  art  of  abstract 
interpretation  is  developing  new  abstract  domains  that  are  tailored  to  a  special  class  of  programs. 
Often  a  new  domain  is  motivated  by  a  new  kind  of  invariant  that  is  need  to  effectively  analyze 
and  prove  properties  of  a  special  class  of  program.  In  this  project,  with  a  special  focus  on  linear 
digital  filters  and  linear  accumulators,  it  is  natural  to  explore  new  abstract  domains  that  generate 
appropriate  invariants  for  those  classes. 


At  the  outset  of  this  project,  our  working  hypothesis  was  that  the  standard  abstract  domains  used 
in  abstract  interpretation  (including  numeric  intervals)  would  give  results  that  are  too  imprecise 
to  be  useful.  Consequently,  the  project  goal  was  to  explore  abstract  domains  that  are  specialized 
to  linear  digital  filters.  As  discussed  later  in  Section  4,  our  working  hypothesis  was  borne  out, 
but,  surprisingly,  as  we  developed  and  applied  more  specialized  domains,  we  discovered 
techniques  that  produce  the  desired  analytic  results  without  the  need  for  abstract  interpretation. 
We  discovered  techniques  for  generating  closed-form  fonnulas  for  LDFs  and  accumulator  codes 
that  yield  sound  and  reasonably  precise  bounds  on  the  output  stream  and  its  floating  point  errors. 
That  is,  we  were  able  to  obtain  our  desired  analytic  results  in  negligible  time. 

3.3.  Integration  of  Model  Analysis  and  Code  Analysis 

This  section  describes  the  original,  standalone  platforms  of  the  two  tools,  and  the  modifications 
and  bridges  that  we  designed  to  allow  their  SANA-enhanced  versions  to  work  as  one. 

3.3.1.  Original  HiLiTE  Platform 

Figure  7  is  an  overview  of  the  HiLiTE  tool  and  its  usage  context,  as  deployed  in  the  avionics 
certification  programs  in  Honeywell  aerospace  products.  HiLiTE  provides  comprehensive  design 
model  analysis  and  test  generation  for  MATLAB  Simulink/Stateflow  models.  It  performs 
symbolic  analyses  combining  computation  semantics,  control  transforms,  and  temporal 
properties  in  a  unified  analysis  framework.  Error  propagation  is  a  recent  addition  to  support 
increasing  tolerances  only  where  necessary  and  identify  when  a  construct  is  untestable  because 
cumulative  error  makes  an  execution  path  undecidable. 
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Significant  reduction  in  certification  cost  &  cycle  time  for  DO-178B/C 


Figure  7.  HiLiTE  Overview 

Range  propagation:  During  model  analysis,  HiLiTE  propagates  signal  dimensions,  data  type, 
and  range  information  from  the  model  input  blocks  through  the  rest  of  the  blocks  in  the  model  to 
the  model  outputs.  In  this  propagation,  the  dimensions,  data  type,  operating  range,  and  hard 
range  constraint  are  determined  for  all  intermediate  signals  (all  inputs  and  outputs  of  each  block 
instance)  in  the  model.  The  range  computation  takes  into  account  the  specific  mathematical  and 
functional  effect  of  each  library  block.  The  following  are  computed  for  each  output  of  each 
block: 

•  Shape:  The  dimensions  of  the  signal  represented  as  a  list  of  sizes  for  each  dimension; 
scalar  signals  have  shape=l. 

•  Data  type:  The  data  type  is  inferred  from  the  block  type  and  input  data  types. 

•  Normal  operating  range:  The  range  of  possible  values  the  output  can  have,  for  all 
combinations  of  input  values  of  the  block  within  the  normal  operating  range  of  each 
input.  This  may  include  an  estimated  error  component  on  each  end  of  the  range. 

•  Maximum  allowable  range  constraint  (also  called  a  “hard  bound”  or  “constraint”):  This 
denotes  a  constraint  on  the  feasible  values  due  to  specific  block  computations.  Examples 
are  outputs  of  the  constant  and  range  limiter  blocks;  values  on  these  signals  outside  the 
constraint  are  not  feasible  in  the  model. 
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Based  on  the  range  analysis  certain  model  design  defects  can  be  analyzed  including  overflow 
conditions,  frozen  signals  (constant  value),  un-testable  conditions,  and  violation  of  modeling 
design  guidelines  or  block-specific  constraints. 

Feedback  loops  and  patterns:  HiLiTE  detects  and  breaks  any  feedback  loops  and  perfonns  the 
propagations  (Shape/DataType/Range)  for  all  the  blocks  interconnected  in  the  loop.  Similarly,  it 
analyzes  certain  combinations  of  blocks  to  identify  common  patterns.  HiLiTE  substitutes  a  single 
block  for  the  model  blocks  that  comprise  the  pattern  and  performs  the  comprehensive 
propagation  analysis  and  test  vector  generation  for  the  entire  pattern  as  a  whole. 

Relationship  propagation:  HiLiTE  has  the  capability  to  detennine  the  relationship  of  a  block’s 
output  to  the  upstream  source  ports  in  the  model.  A  constraint  in  the  model  occurs  when  a  signal 
fans  out,  creating  multiple  branches  that  converge  downstream  onto  the  inputs  of  a  single  block. 
When  two  or  more  such  inputs  are  totally  related  to  each  other  (identity  or  linear  relationship 
with  zero  slope)  HiLiTE  recognizes  this  and  maintains  a  list  of  totally  related  inputs  for  the 
blocks  in  the  model.  If  the  output  of  the  block  is  a  polynomial  of  the  block  inputs  HiLiTE 
recognizes  the  polynomial  relationship  and  uses  it  to  determine  a  tight  range  bound  at  the  output. 

The  relationship  propagation  takes  into  account  the  specific  mathematical  and  functional  effect 
of  each  library  block.  The  propagation  of  these  relationships  supports  the  analysis  of  the  data 
flow  expression  to  recognize  and  solve  the  equation  to  derive  tight,  exact  range  bounds  at  a 
block’s  output.  Also  it  helps  in  identifying  test  generation  constraints  due  to  related  inputs  in  the 
model. 

Feedback  loop  invariant  analysis:  We  extended  HiLiTE’s  feedback  loop  and  relationship 
processing  to  discover  cycle  invariants  describing  cycle  input-state-output  behaviors.  A  cycle 
invariant  is  expressed  by  guard/assignment  pairs,  where  each  pair  represents  a  possible  execution 
path  of  the  cycle.  HiLiTE  discovers  cycle  invariants  by  performing  path  analysis  and  integration 
of  individual  block  invariants  following  the  analysis  procedures  described  below. 

Figure  8  shows  a  model  of  a  limited  counter  with  only  one  cycle.  This  cycle  has  three  possible 
paths  due  to  the  three  different  guard  conditions  of  block  limit.  For  each  path,  we  collect  the  set 
of  guard/assignment  pairs  from  blocks  within  the  cycle  starting  with  the  time-dependent  block 
delay,  followed  by  sum  and  limit. 


■■I. - - 

UpCounlerLimnSlop  (i)  [1]  | 

Pathl:  counter  input  lower  than  lower  bound. 

Guard:  y(k)_sum  <  0,  Assignment:  y(k)_limit  =  0 
Path2:  counter  input  between  lower  and  upper  bound. 

Guard:  0  <=  y(k)_sum  <  =  9,  Assignment:  y(k)_limit  =  y(k)_sum 
Path3:  counter  input  greater  than  upper  bound. 

Guard:  y(k)_sum  >  9,  Assignment:  y(k)_limit  =  9 
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Figure  8.  Simulink  model  of  limited  counter  and  three  possible  paths 

Next,  we  use  SMT  solver  Z3  to  check  the  satisfiability  of  each  path.  HiLiTE  automatically 
generates  the  input  file  for  Z3  containing  variable  declarations,  sets  of  constraints  translated  from 
the  block  level  guard/assignment  pairs  and  some  checking  commands.  Z3  returns  “unsat”  for 


Approved  for  Public  Release;  Distribution  Unlimited. 
16 


each  invalid  path,  which  is  then  removed  from  the  cycle  invariant.  For  the  remaining  valid  paths 
(all  3  paths  are  valid  in  Figure  8),  we  further  reduce  the  block  level  guard/assignment  pairs  into  a 
cycle  level  guard/assignment  pair  possessing  only  the  cycle  level  input/output  variables.  This 
process  includes  virtual  substitution  to  eliminate  the  intermediate  variables  as  well  as  some 
general  simplification.  Finally,  the  cycle  invariant  is  printed  out  in  the  log. 

The  OptionSet  information  to  include  in  the  command  file  looks  like  this: 

<OptionSet  name="ModelAnalysisDirectives"> 

<Option  key="  AnalyzeFeedbackLoops ">True</Option> 

</OptionSet> 


By  running  HiLiTE  on  a  model  with  the  option  "AnalyzeFeedbackLoops"  set  as  “True,”  HiLiTE 
obtains  the  cycle  invariant  and  decoupler  outport  range.  In  the  analysis  process,  messages  are 
given  as  information  including  “cycle  list  information,”  “path  expressions,”  and  “range 
calculated  for  the  time-dependent  block  outport  in  current  cycle  list”. 

3.3.2.  Original  CodeHawk  Platform 

Kestrel  Technology’s  CodeHawk  technology  is  an  Abstract  Interpretation  framework,  written  in 
Ocaml,  that  is  meant  to  be  specialized  to  particular  programming  languages  and  program 
properties  of  interest.  A  specific  analyzer  tool  is  built  for  a  specific  combination  of 
language/properties  by  specializing  the  Ocaml  source  and  compiling  to  the  intended  operating 
system  (Figure  9  below).  Thus  CodeHawk  can  in  principle  run  on  any  machine/OS  for  which 
there  is  an  Ocaml  compiler.  We  have  standard  front  ends  for  C,  Java  and  x86  binaries,  and  a 
variety  of  abstract  domains  over  which  to  approximate  program  behaviors.  For  this  project,  we 
started  with  an  exiting  specialization  for  memory  safety  properties  of  C  programs  that  typically 
runs  on  Unix/Linux/MacOS.  The  need  to  extend  this  analyzer  to  a  new  floating-point  abstract 
domain  during  the  course  of  this  project  caused  us  to  re-implement  a  portion  of  the  analyzer  in 
C++,  thus  breaking  the  monolithic  Ocaml  compilation  and  resulting  in  a  mixed  language  system 
that  (currently)  runs  only  on  MacOS  (see  3.3.4  below). 
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3.3.3.  JSON  Interchange  Architecture 

Since  HiLiTE  and  CodeHawk  typically  run  on  different  operating  system  platforms,  we  made  an 
early  architectural  decision  to  exchange  infonnation  between  their  SANA-enhanced  versions  via 
an  external  file  interface.  This  would  allow  development  to  proceed  in  parallel  at  different  sites, 
and  for  the  final,  integrated  deliverable  to  be  run  on  the  same  or  different  machines  at  the  same 
site,  as  well  as  different  machines  at  different  sites.  We  designed  a  custom  JSON  interface 
specification  that  was  sufficiently  abstract  to  mask  differences  between  the  two  tools  and  their 
respective  source  languages.  As  the  first  year  enhancements  and  integration  progressed,  we 
repeatedly  revised  this  common  JSON  specification  to  match,  buffering  the  changes  with 
grammars  and  parser  generators  to  rebuild  the  necessary  translators. 

Figure  10  shows  the  overall  file  exchange  integration  for  Linear  Digital  Filters.  CodeHawk 
normally  takes  as  input  all  of  the  C  sources  required  to  compile  the  filter  application  (or  any 
other  C  application).  For  this  project,  we  restricted  the  C  source  to  just  that  implementing  the 
actual  filter  function  that  HiLiTE  generates.  In  addition  to  the  C  sources,  HiLiTE  also  passes 
model-level  information  about  the  filter  (parameters,  recurrence  relation,  input  bounds,  etc)  in  a 
JSON  file,  which  contains  a  reference  to  the  C  source  file.  CodeHawk  consumes  both  inputs, 
uses  the  model-level  JSON  information  to  specialize  the  abstract  interpretation  to  the  input  filter, 
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then  writes  back  an  approximation  of  the  value  and  error  bounds  of  the  filter’s  output  to  the  same 
JSON  instance.  HiLiTE  then  consumes  this  rewritten  instance  to  exploit  the  computed  bounds  at 
the  filter  model  level. 


Integrated  Model  Level  and  Code  Level  Analysis 


for  each  filter  instance.  bounds 


3 


Figure  10  HiLiTE  <->  CodeHawk  Integration  via  files 
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File  Exchanged  forbetween  Tools 

{  "filterlnstances" :  [ 

{"objectType":  "FilterlnstanceSpec", 

"instancelD":  1, 

"filterType":  "ResetLeadLagFilterFixedTau", 

"cFile":  "RefMdIResetLeadLagFilter.c", 

"initFunction":  "RefMdlResetLeadLagFilterlnit", 

"computeFunction":  "RefMdl_ResetLeadLagFilter_p", 

"inputs"  :  [{"u":  "rtu_0”}], 

"outputs":  [{"y":  "rty_0"}], 

"paramAssoc":  [{"rate":  "rtp_rate”},  {"Td":  "rtpTd"},  {"Tn":  "rtpTn"}], 
"paramValues":  [{"rate":  40),  {"Td":  1),  {"Tn":  1}], 

"inputRanges":  [  {"u" :  {"objectType" :  "Range",  "type":  "float", 

"intervals":  [{"min":  {"inclusion":  true,  "value":  -100}, 
"max":  {"inclusion":  true,  "value":  100  }})} 

}, 

{"r" :  {"objectType" :  "Range",  "type":  "bool", 

"intervals":  [{"min":  {"inclusion":  true,  "value":  0}, 

"max":  {"inclusion":  true,  "value":  0  }}]} 

}, 

{"rv" :  {"objectType" :  "Range",  "type":  "float", 

"intervals":  [{"min":  {"inclusion":  true,  "value":  0}, 

"max":  {"inclusion":  true,  "value":  0  }}]} 

}  1. 

"outputRanges":  [  {"y"  :  {"objectType" :  "Range",  "type":  "float", 

"intervals":  [{"min":  {"inclusion":  true,  "value":  -200}, 
"max":  {"inclusion":  true,  "value":  200}})} 

}  ], 

"errorBoundFuncParams":  {"m" :  1,  "b" :  .01} 

}, 

{"objectType":  "FilterlnstanceSpec", 

"instancelD":  2,  26 


Figure  11  HiLiTE  <->  CodeHawk  JSON  Exchange  Example 

Figure  1 1  illustrates  a  sample  JSON  input  file  (portion  thereof)  for  a  filter  analysis.  This  file 
contains  information  about  the  filters  found  in  a  model.  A  given  C  application  may  contain  a 
number  of  filter  functions,  and  any  one  filter  function  may  be  used  by  multiple  “instances”  of 
that  filter.  The  JSON  file  for  filter  analysis  consequently  uses  a  series  (JSON  array)  of 
filterlnstanceSpec  objects.  Each  instance  represents  a  particular  set  of  parameters  for  the  filter, 
theoretical  bounds  on  the  inputs  to  the  filter,  references  to  the  C  sources  and  their  algorithmic 
parts  (init  and  step  functions)  that  implement  the  filter,  and  placeholder  JSON  objects 
representing  the  value  and  error  bounds  on  the  outputs  from  the  filter.  CodeHawk  will  consume 
the  input  characteristics,  and  rewrite  the  JSON  output  objects  to  record  the  results  of  analysis. 

Note  that  the  HiLiTE  JSON  output  refers  to  other  files.  These  files  may  be  the  source  code  for 
standard  functionality  referenced  by  the  analyzed  model  or  it  may  be  the  source  code  generated 
for  the  model  itself. 

The  JSON  format  for  integrating  accumulators  explored  in  the  second  year  of  this  project  is 
different  than  that  for  filters  since  different  infonnation  is  required  for  these  two  types  of 
analyses.  We  avoid  confusion  in  the  exchange  both  by  using  different  files  for  each  type  of 
analysis  data  and  by  the  creation  of  a  ‘filterlnstances’  list  for  filter  analysis  data  and  an 
‘accumulatorlnstances’  list  for  accumulator  analysis  data. 
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3.3.4.  Tool  Software  Integration  Architecture 

The  toolset  integration  has  several  ‘levels’.  There  is  the  integration  of  the  entire  system  of  tools 
(CodeHawk  +  HiLiTE)  as  well  as  the  integration  of  each  of  the  tool’s  internal  systems  within 
itself. 

HiLiTE  is  a  Windows  .NET  application  comprised  of  multiple  C#  and  J#  assemblies,  some 
statically  linked  and  others  dynamically  loaded.  In  addition,  HiLiTE  links  to  COM  libraries  for 
its  interface  to  MATLAB  Simulink  to  open  and  process  the  contents  of  a  model. 

CodeHawk  is  normally  an  OCaml  application  specialized  to  a  particular  source  language  (C, 
Java,  x86  binary)  and  a  particular  set  of  program  conjectures  to  be  proved.  This  project  required 
specializing  Kestrel’s  existing  C  analyzer  to  a  new  abstract  domain  of  floating  point  numbers. 
The  libraries  for  this  domain,  however,  are  written  in  C++,  so  our  original  architectural  plan  was 
to  write  Ocaml  “wrapper”  code  for  these  C++  libraries  to  bring  them  into  the  existing  Ocaml 
architecture.  This  proved  to  be  infeasible,  so  we  were  forced  to  write  an  entirely  new  abstract 
interpretation  engine  in  C++  and  split  the  erstwhile  monolithic  CodeHawk  architecture  into 
several  pipelined  pieces.  Because  of  the  exigencies  of  developing  and  integrating  these  pieces 
from  multiple  languages  at  multiple  sites,  the  resulting  SANA-specialized  CodeHawk  (currently) 
only  runs  reliably  under  MacOS.  Its  constituent  pieces  are  assembled  from  C,  C++,  and  Ocaml 
compiled  binaries,  plus  a  Unix  shell  script,  and  integrated  into  a  single  application  as  a  runnable 
Java  jar  file. 

A  consequence  of  this  mixed-language/mixed-operating-system  configuration  is  that  the 
combined  HiLiTE/CodeHawk  system  typically  runs  on  separate  Windows  and  Mac  machines, 
with  round-trip  communication  implemented  via  shared  C  source  and  JSON  files. 

The  source  code  for  the  library  of  standard  functionality  that  can  be  referenced  is  installed  as  part 
of  the  CodeHawk  installer.  Currently  this  set  of  standard  functionality  consists  of  all  the  known 
filters. 
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Figure  12  HiLiTE  <->  CodeHawk  Comm  Process 

The  system  of  tools  is  integrated  through  the  use  of  the  files  described  above.  Since  these  JSON 
files  also  reference  ‘C’  source  files,  those  files  must  either  already  be  present  on  the  CodeHawk 
machine  or  must  be  transferred  along  with  the  JSON  file  that  references  them.  The  filter  source 
code  is  installed  along  with  CodeHawk  and  thus  need  not  be  transferred.  HiLiTE  recognizes  this 
and  bundles  model  specific  source  code  along  with  the  JSON  in  a  package  that  can  be  easily 
made  accessible  to  the  CodeHawk  machine. 

Because  the  HiLiTE  machine  already  has  all  the  source  code,  the  only  information  that  need  be 
transmitted  back  from  the  CodeHawk  machine  is  the  result  JSONfile. 

The  actual  mechanism  for  transfer  can  be  through  use  of  a  file  server,  email  or  any  other 
mechanism  convenient  to  the  user.  As  part  of  moving  files,  they  will  need  to  be  placed  where 
CodeHawk/  HiLiTE  can  find  them.  See  the  User  Guide  for  details. 
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4.  RESULTS  AND  DISCUSSION 


4.1.  Discoveries  in  Abstract  Interpretation 

As  discussed  in  the  previous  section,  our  approach  was  to  explore  various  abstract  domains 
within  an  abstract  interpretation  framework,  implemented  by  CodeHawk,  for  analyzing  linear 
digital  filters  and  accumulating  integrators.  The  analytic  results  that  we  desire  are 

(1)  bounds  on  the  output  of  a  filter, 

(2)  bounds  on  the  floating-point  errors  in  the  output  of  a  filter, 

(3)  bounds  on  the  output  of  an  accumulating  integrator  as  a  function  of  the  number  of 
iterations,  and 

(4)  bounds  on  the  floating-point  errors  in  the  output  of  an  accumulating  integrator  as  a 
function  of  the  number  of  iterations. 

4.1.1.  Interval  Abstract  Domain 

As  always  in  analyzing  a  program,  the  first  question  is  what  kinds  of  invariants  are  required  to 
establish  the  desired  results.  Since  we  obviously  want  bounds  as  results,  the  first  abstract 
domain  to  try  is  intervals. 

As  a  running  example,  we  will  use  a  simple  first-order  linear  digital  filter  (called  a  Lead  Lag 
filter,  see  Section  4.2.1): 


9  7  11 

Yn+l  ~  ^2  ^n+1  —  13  Un  13  yn 


(3) 


where  the  input  stream  is  written  u0,  ulr ...  and  the  output  stream  is  y0,  ylt ...  We  are  given  that 
the  range  of  the  input  values  is  [-1,1],  and  that  y0  =  0.  What  are  the  possible  values  for  the 
output  stream?  If  we  use  the  interval  abstract  domain,  representing  y0  as  [0,0]  and  iterate  the 
recurrence  we  get 


9  7  11 

Vi  =  777  [-U]  -  —  [-1,1]  +  —[0,0] 
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13 


13 


=  [- 


-16  16 


13  '13 


3/2  13  ^  1'1^  13  ^  + 


11 

(-16  16] 

(-192  192] 

13 

[  13  '  13J 

[  169  '169.1 

(4) 

(5) 


Note  that  the  sequence  is  growing  and  we  have  gone  through  the  iterative  loop  twice.  A  typical 
step  in  abstract  interpretation  is  to  speed  up  convergence  to  a  fixpoint  by  applying  a  widening 
operator.  There  are  standard  widening  operators  for  each  abstract  domain  and  their  choice  is 
somewhat  of  an  art.  A  widening  operator  suggested  in  earlier  work  on  filters  [3]  is  to  jump  up 
to  the  nearest  power  of  ten,  so  [-10, 10]  in  our  case.  This  does  in  fact  converge  since 

9  7  11 

—  [—1,1]  -  —[—1,1]  +  —  [—10,10]  c  [-10,10],  (6) 
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The  actual  range  of  the  output  is  [-1,1],  so  this  inferred  range  is  very  imprecise,  although  sound. 
The  best  that  can  be  obtained  using  the  interval  domain  is  [-8,8]  which  is  better,  but  far  from  the 
exact  bound. 


4.1.2.  Cancellation  Abstract  Domain 

The  essential  reason  for  the  imprecision  when  using  the  interval  domain  is  that  there  are 
cancellation  possibilities  if  we  simply  unroll  the  recurrence: 


9  7  11 

y2=T3“2“l3U‘+T3yi 


11  r  9 


=  - u-1  —  - Ui  +  — 

13  2  13  1  13 


11 


- Mi  —  - Un  +  - Vn 

13  1  13  0  13  "Xo 


—  T7rM.2  +  Mi 

13  2  169  1 


77 

169 


M0 


(V) 


where  some  of  the  contribution  of  ul  to  y2  has  been  cancelled  out.  This  phenomenon  suggests 
an  approach  to  constructing  an  abstract  domain  that  generates  stronger  invariants  for  linear 
recurrences:  supplement  intervals  with  symbolic  expressions  to  allow  cancellation  as  part  of  the 
abstract  operators  of  the  domain.  As  an  example,  instead  of  abstracting  y2  to  an  interval,  we 
abstract  it  to  a  symbolic  sum  (where  Abs(x)  denotes  the  abstract  value  of  variable  x ): 


Abs(y2)  =  ^ Abs(u2 )  +  [-^9 ’1^9]  (8) 

which  will  allow  us  to  partially  cancel  Abs(u2 )  in  the  next  iteration  (to  generate  Abs(y?i)).  We 
implemented  an  abstract  domain  of  this  type  and  indeed  it  converges  to  the  exact  bound  [-1,1]  in 
our  example.  It  generally  performs  well  for  first-order  linear  filters  since  there  are  a  fixed 
number  of  terms  to  be  carried  between  iterations  to  harvest  all  the  cancellation  possibilities. 
Unfortunately,  for  higher-order  linear  filters,  there  is  no  fixed  bound  on  the  number  of  terms  to 
be  carried  so  this  approach  won’t  work  as  stated.  We  explored  several  ways  to  further  elaborate 
the  approach  for  second-order  filters,  but  the  result  is  very  complex  and  of  unclear  benefit. 


Fortunately,  the  idea  of  unrolling  the  recurrence  carries  the  seed  of  an  alternate  approach  that 
proved  successful.  Intuitively,  we  can  imagine  fully  unrolling  the  recurrence  until  any  output  yn 
is  expressed  as  a  linear  sum  of  all  preceding  inputs  (i.e.  all  previous  outputs  in  the  recurrence  are 
eliminated  by  unrolling).  Before  we  turn  to  the  approach  based  on  that  concept,  we  first  review 
our  parallel  efforts  to  analyze  the  error  due  to  floating-point  representations. 


4.1.3.  Analyzing  for  Error  Bounds 

In  addition  to  wanting  bounds  on  the  values  of  the  output  of  a  filter  or  accumulator,  we  also  want 
to  know  bounds  on  the  accumulated  errors  due  to  floating  point  implementation  of  real  values. 
Using  simple  intervals  for  error  analysis  will  result  in  a  non-converging  series  of  bounds.  In 
practice  the  errors  are  bounded  and  fairly  small,  at  least  for  first-order  filters,  so  once  again  the 
interval  domain  proves  to  be  too  imprecise  for  this  application. 
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As  mentioned  earlier,  if  the  standard  abstract  domains  do  not  provide  a  useful  time/precision 
tradeoff,  then  the  skilled  abstract  interpretationist  explores  the  type  of  invariants  that  are  required 
to  obtain  the  desired  results  and  develops  abstract  domains  to  generate  those  invariants. 

We  define  the  error  as 


f(x,ri)  —  r(x,n )  <  e(x,ri)  (9) 

where  r(x,  n)  denotes  the  real  value  of  variable  x  at  the  nth  iteration,  /  (x,  n)  is  the  floating  point 
value  of  x  after  n  iterations,  and  e(x,  ri)  is  the  error,  expressed  as  a  bound  on  the  difference  of  the 
two. 

If  we  view  the  output  yn  as  a  linear  sum  of  all  preceding  inputs,  we  get  a  sense  that  the  net  error 
due  to  floating  point  can  be  calculated  as  a  linearly  weighted  sum  of  the  error  contributions  of 
each  of  the  inputs  and  the  operations  performed  upon  them.  This  leads,  after  some  trial  and 
error,  to  seeking  an  invariant  of  the  form 

e(x,n)  =  Ax  +  MxbPxYn  ?  (10) 

^~Ji=0 


where  e(x,  ri)  is  the  error  range  of  variable  x  at  a  code  location.  The  sum  conveys  the  effect  of 
the  errors  from  all  earlier  inputs.  The  invariant  is  parametric  on  A,  M,  p,  b,  and  /?,  where  Ax,  Mx, 
and  px  vary  per  variable  x.  As  we  simulate  the  concrete  program  (implementing  the  recurrence) 
via  abstract  operators,  the  effect  will  be  to  update  the  current  abstract  values  of  Ax,  Mx,  and  px  . 

We  have  been  able  to  derive  the  abstract  operations  that  preserve  the  form  of  the  invariant  by 
changing  the  parameters.  To  illustrate,  consider  the  case  that  the  concrete  program  performs  a 
sum  of  two  variables 


v  —  w  +  z 


(ID 


where  the  invariant  holds  before  the  assignment.  We  want  to  calculate  updates  to  the  parameters 
of  the  invariant  such  that  the  invariant  holds  after  the  assignment. 

We  assume  that 


fr(w)  =  Aw  + 


e(z)  =  Az  + 


Zn 

Pi 

i= 0 


M  UPw 

‘v‘wu  /  rw 


MzbPz 


Zn 

Pi 

i= 0 


(12) 


(13) 


then  (letting  em  denote  the  machine  error  for  the  chosen  floating-point  representation) 
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f(v,  n)  —  r(v,  ri) 

—  (/(w, n)  +/(z,n))(  1  +  8)  —  (r(w,n)  +  r(z,n))  where  |5|  <  em 

—  (r(w,  ri)  +  e(w)  +  r(z,n)  +  e(z))(l  +  8)  —  (r(w,n)  +  r(z,n)) 

=  (r(w,n)  +  r(z,n))5  +  (e(w)  +  e(z))(l  +  5) 

=  r(v,n)<5  +  (e(w)  +  e(z))(l  +  <5)  (14) 

<  max(r(u,  n))  em  +  (Aw  +  MwbPw  Yi=o  PI  +  +  MzbPz  If=0  ^)(1  +  8 ) 

<  [max(r(u,  n))  +  (i4w  +  i4z)(l  +  em)] 


+ 


(Mw  +  Mz)(  1  +  fm) 
b 


£,max(pw,pz)+l 


=  e(v,n) 


where  we  update  the  parameters  according  to 

Av  =  [max(r(v,n))  em  +  ( Aw  +  Az)(l  +  em )]  (15) 

=  (Mh,  +  Mz)(  1  +  em) 
v  b 

pv  =  max(pw,pz)  +  1.  (17) 

That  is,  when  the  concrete  implementation  assigns  a  sum  of  two  variables  to  a  variable,  we 
perfonn  the  updates  above  in  the  abstract  domain.  Similar  abstract  operators  have  been  derived 
and  implemented  for  other  concrete  operations:  subtraction,  multiply  by  a  constant,  divide,  and 
so  on. 


Using  this  abstract  domain  resulted  in  finite  error  bounds  that  were  surprisingly  small.  One 
might  intuit  that  floating-point  errors  should  accumulate  without  bound  as  a  linear  recursion 
progresses.  However,  it  seems  that  filters  inherently  have  a  discounting  mechanism  so  that  the 
contribution  of  previous  inputs  increasingly  diminish  with  age.  The  same  mechanism  seems  to 
apply  to  errors:  the  effects  of  older  errors  are  discounted  with  age,  so  that  the  age -weighted 
errors  form  a  series  that  converges. 

4.1.4.  Closed-form  Solutions 

The  explorations  discussed  in  previous  subsections  pointed  to  the  conceptual  value  of  fully 
unrolling  the  recursion  and  expressing  the  current  output  yn  as  a  linear  sum  of  all  preceding 
inputs.  Anca  Browne  discovered  that  the  coefficients  of  that  linear  sum  are  fixed  by  the  linear 
recurrence,  and,  moreover,  are  themselves  characterized  by  a  homogeneous  linear  recurrence 
relation.  This  remarkable  fact  allows  us  to  find  a  closed-form  expression  of  the  coefficients  and 
then  to  symbolically  calculate  exact  bounds  on  the  output  values.  A  generalization  of  this 
approach  can  also  be  applied  to  calculating  bounds  on  error  values  and  to  find  bounding 


Approved  for  Public  Release;  Distribution  Unlimited. 
26 


functions  on  the  range  and  errors  of  accumulating  integrators.  The  effect  is  that  through 
symbolic  analysis  of  the  linear  recurrences  that  characterize  linear  digital  filters  and 
accumulating  integrators,  we  can  derive  exact  or  tight  bounds  on  the  range  of  outputs  and  their 
floating  point  errors,  in  essentially  constant  time.  This  stands  in  contrast,  say,  to  20  hours  of 
analysis  time  needed  to  perform  sophisticated  abstract  interpretation-based  analysis  for  just  one 
Airbus-related  analysis  in  [3]. 

Let 

Tn+1  =  aOun+l  +  alun  +  ■  ■  ■  +  aNun-N+l  +  ^lTn  +  ■  ■  ■  +  ^wTn-W+1  O8) 

be  the  recurrence  equation  for  the  output  of  a  digital  filter  of  order  N  that  defines  the  output  yn+1 
in  terms  of  N  +  1  present  and  past  inputs  and  N  past  outputs.  Unfolding  the  equation,  yn+1can 
be  expressed  as  a  linear  combination  of  un+1  ...  u^. 

y-n+1  =  v+i  +  c?+ v  + ...  +  crv  (i9) 

Notice  that  for  any  n,  cj}+1  —  a0.  Therefore  for  any  n,  cf+1=  ax  +  hxCo  =  ax  +  b±a 0,  etc. 

The  key  point  to  notice  is  that  c£  does  not  depend  on  n,  for  any  k.  Without  the  superscript, 

Tn+l  ^0^-n+l  +  Ci^-n  T  •  •  •  T  (20) 

The  first  N+l  terms  c0, ... ,  cN  depend  on  both  a0, ... ,  aN  and  blt ... ,  bN.  However,  for  any  n  >  N, 

Ui+l  =  Vri  +  ^2Cn-l  +  ■■■+bNcn_N+1  (21) 

This  homogeneous  linear  recursion  fonnula  of  order  N  has  the  general  solution 

Ui+i  =  Pi(.ri)ri  +  P2(n)r 2n  +  ...  +  P;(n)ryn  (22) 

where  r1( ... ,  ry  are  solutions  of  the  characteristic  equation 

xN  =  b1xN~1  -I-  b2xN~2  +  ...  +  bN  (23) 

with  respective  multiplicities  m v  m2, . . . ,  and  Pj's  are  polynomials  of  degree  mi  —  1.  See 

e.g.  [4]  for  more  background  on  solving  (in)homogeneous  linear  recurrence  equations. 

The  initial  conditions  that  detennine  the  solution  are  given  by  cv  ... ,  cN. 


4.1.5.  Range  Analysis  for  First-order  Linear  Digital  Filters 

The  general  fonn  of  a  first-order  linear  filter  is 
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y0  —  o 

yn+ 1  =  a0Wn+i  +  aiun  +  bryn  for  n  >  0 

where  it0  =  0.  The  coefficients  ct  are  given  by 

c0  —  ao 

ci  —  ai  +  Mo 

cn+i  =  C\K 

and  therefore  the  output  is 

yn+ 1  =  a0un+1  +  (%  +  MoH^n  +  bx  un_!  +  •••  +  MM) 


(24) 

(25) 

(26) 

(27) 

(28) 

(29) 


or 


n-l 


yn+ 1  =  a0un+1  +  (%  +  Mo)  ^  Mn-i 

i= 0 


(30) 


Since  we  are  assuming  that  all  inputs  have  the  same  range,  we  can  omit  the  subscript  on  inputs 
for  purposes  of  computing  the  range  bounds  of  the  output: 


n-l 


max(yn+1 )  =  max(a0it)  +  max((a1  +  bxa 0)  ^  b[u) 

i= 0 


1  -b\ 


—  max(u)(a0  +  ^  (ax  +  b-^ao)  ■  i_  ^  j) 

1 


<  max(it)  a0  +  (ax  +  /^cio) 


1  - 


(31) 


when  <  1.  Similarly, 


min(yn+ 1)  <  min(u )  a0  +  (ax  +  f^cio) 


1  ~  M 


As  an  example,  consider  the  lead  lag  filter  (see  Section  4.2.1)  described  by 

To  —  0 

9  7  11 

yn+l  —  ^2  ^ n+1  —  ^  g  ^  3  ~^n 


(32) 

(33) 
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with  u0  —  0  and  for  any  i  >  0,  iq  6  [—1, 1].  Then  the  coefficients  c*  are  given  by 


Ci  - 


+ 


9 

13 

11  9  8 


13  13  13  169 

8  11  n 
Cn+i  ~  ( fj) 


and  the  range  of  the  output  by 


9  8 

max(yn )  <  — —  + 


13  '  169-11 
1  13 


min(yn)  <  -  — 
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13  169  -  _  11 

1  13 


=  1 


=  -1 


(34) 

(35) 

(36) 


(37) 

(38) 


4.1.6.  Range  Analysis  for  Second-order  Linear  Digital  Filters 

A  similar,  slightly  more  complex  treatment  computes  the  range  bound  for  second-order  linear 
filters: 


y  0  =  yi  =  0 

yn+1  =  a0un+1  +  atun  +  a2un_1  +  hxyn  +  h2yn_1  for  n  >  1. 
Suppose  that  the  characteristic  equation  has  two  distinct  solutions  rt  and  r2 .  Then 

c0  —  ao 

C-y  —  CLy  “f  t)yCLQ 
c2  =  by  +  bfa0  +  b2a0 
Cn+l  =  sri  +  tr2n 

where  s  and  t  are  the  solutions  of 

Cy  —  s  +  t 

c2  =  sry  +  tr2 


(39) 

(40) 


(41) 

(42) 

(43) 

(44) 

(45) 

(46) 
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The  output  is 


yn+ 1  =  a0un+ 1  +  (s  +  t)un  +  (srj  +  tr2)un_1+...  (47) 


In  the  case  when  s,  t,  rlt  r2  are  all  positive  the  range  is  given  by 


max(yn )  —  max(a0u)  +  s(max(u))(l  +  +  r2  +  . . . )  +  t(max(it))(l  +  r2 

+  r|  +  . . . ) 

min(yn )  =  min(a0  u)  +  s(mm(it))(l  +  +  r12+...)  +  t(min(it))(l  +  r2 

+  '•!  +  ...) 


(48) 

(49) 


For  example,  consider  the  following  quadratic  fdter  (discussed  in  Section  4.2.4) 

To  =  y±  =  0  (50) 

592841  i  1185682  _  592841  _  126814318  51175081 

yn+1  _  78010601  Un+1  +  78010601  Un  +  78010601  Un_1  +  78010601  yn  ~  78010601  yn_1 

(51) 

where  u0  —  0  and  for  any  i  >  0,  ut  E  [—100, 100], 

The  first  coefficients  are 


592841 

cO  -  - 

(52) 

78010601 

167676492512320 
cl  -  - 

(53) 

6085653868381201 

22504866023935154502080 

(54) 

474745515750392387111801 

and  the  others  are  determined  by  the  solutions  of  the  characteristic  equation 

2  _  126814318  51175081 

X  “  78010601  X  ~  78010601 


(55) 


which  are 
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63407159  +  160  VH04257321 
r12  =  -  *  {0.744646,  0.880957} 


(56) 


The  general  solution  is 


(57) 


-n+1 


=  s  (■ 


63407159  -  160 V110425732T  63407159  +  160  a/11042573211 

- )n  +  tt - )n 

'  v  700 1  r\  J 


78010601 


7801060 


where  s  and  t  are  solutions  of  the  system 


s  +  t  =  c-. 


s  (- 


63407159  -  160a/1104257321  63407159  +  160  a/1104257321 

- )  +  t  (• - )  =  c7 

'  v  7001  o/:oi  J  z 


78010601 


78010601 


(58) 

(59) 


which  are 


(60) 


s 


t 


5579860  (—976009979996381  +  391157262321  VH04257321) 
39530163748423009547191 

5579860  ^64956313  (391157262321  +  883861  VH04257321) 

6085653868381201  * 


*  -0.169695 


(61) 

0.197247 


The  expansion  of  yn+1  is 


yn+ 1  =  c0un+1  +  (s  +  t)  itn  +  (srj  +  t  r2)itn_1  +  •••  +  (s  rf  1  +  t  r2n  1)u1  (62) 

Notice  that  s  +  t  r2  >  (s  +  t)r2fc  >  0  since  s  <  0,  s  +  t  >  0  and  r2  >  rx.  As  all  the 
coefficients  are  positive,  the  maximum  is  attained  for  max(it)  and  the  minimum  for  min(u ) 


max  (yn)  =  100  (c0  +  (s  +  t)  +  (s  rx  +  t  r2)  +  . . . ) 

1  1 


—  100  (c0  +  s  - - 1-  t- - )  —  100 

1  —  1  —  r2 


(63) 


min  ( yn )  =  —100  (c0  +  (s  +  t)  +  (s  ?i  +  t  r2)  +  . . . ) 

—  — 100  (^Cq  +  5  — -  +  t  — - 'j  —  — 100. 

V  1  -  ry  1  -  r2J 


(64) 
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4.1.7.  Digital  Filters:  Error  Range 


The  filter  recursive  formula  does  not  immediately  translate  into  a  constant-term  recursive 
fonnula  for  the  floating-point  output  or  its  accumulated  error.  Different  inputs,  with  different 
errors,  will  influence  the  constants  differently.  The  solution  we  adopt  is  to  use  intervals. 

However,  just  reproducing  our  approach  for  finding  the  real  bounds  does  not  work  because  the 
characteristic  equation  cannot  be  solved  in  the  same  way  over  intervals. 

As  before,  the  error  for  variable  x  is 

f{x,  n)  —  r(x,n)  <  e(x,n)  (65) 

where  r(x,  ri)  denotes  the  real  value  of  variable  x  at  the  nth  iteration  for  some  particular 
sequence  of  inputs,  f(x,  n )  is  the  floating  point  value  of  x  after  n  iterations,  and  e(x,  n)  is  the 
error,  expressed  as  a  bound  on  the  difference  of  the  two.  We  will  assume  that  all  the  variables 
are  bounded  over  all  iterations  and  we  will  denote  by  R(x)  an  interval  that  includes  all  possible 
real  values  for  x  and  by  EV  (x)  an  interval  that  includes  all  the  possible  errors  for  variable  x.  For 
inputs  itj,  we  will  take  EV (uL)  to  be  a  symmetric  interval. 

4. 1.7.1.  Symbolic  Interval  Combinations 

Let  J  be  the  set  of  intervals  over  real  numbers.  Let  S  be  a  set  of  symbols  and  <A  be  the  set  of 
affine  combinations 


Gy  Si  +  . . .  Gp  Sp  +  G  (66) 

where  G0, . . .  Gp,  G  G  0  and  S0, . . . ,  Sp  G  S.  Two  combinations  that  differ  only  by  the  order  of 
their  terms  or  by  a  zero-coefficient  tenn  as  in 

Gy  Sy  +  ...  Gp  Sp  +  G  =  Gy  Sy  +  ...  Gp  Sp  +  [0,0]Sp+i  +  G  (67) 


are  considered  identical. 

We  will  assume  that  S  contains  all  the  symbols  that  we  need,  among  which  are  symbols  for  the 
infinite  sequence  of  input  errors  and  output  errors. 

We  define  the  following  operations  on  c/Z: 

1 .  multiplication  by  an  interval 

(GySy  +  ...GpSp  +  G)  *  I  =  (Gy*I)Sy  +  ...(Gp*I)Sp  +  G*I  (68) 

2.  division  by  an  interval 

GySy  +  ...Gp  Sp  +  G)/I  —  (Gy/I)  Sy  +  ...(Gp/I)  Sp  +  G  /  /  (69) 
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3.  addition 


(Gi5i  +  -GpSp  +  G)  +  (G'1S1+  ...G'pSp  +  G')  ^ 

—  (Gi  +  G'i)  +  ■  ■  ■  (Gp  +  G'p)  Sp  +  (G  +  G') 

Notice  that  we  can  always  assume  two  combinations  have  the  same  symbols  because  we 
can  always  add  zero-coefficient  terms. 

4.  subtraction 


{Gy  Si  +  ...  Gp  Sp  +  G)  —  (G'i  Si  +  . .  .G'  p  Sp  +  G 
=  (Gy  —  G’y)  Sy  + ...  (Gp  —  G'p)  Sp  +(G-G') 

5.  Simultaneous  substitution.  Here  is  the  result  of  one  symbol  being  substituted: 

Subst(  Gy  +  ...  Gp  Sp  +  G ;  HySy  +  . . .  HpSp  -I-  H 

—  Gy  (HySy  +  . . .  HpSp  +  +  G2  S2  +  •••  Gp  Sp  +  G 


(71) 


(72) 


Any  V:  S  ->  J  can  be  extended  to  V:  A  ->  J  by 

V(Gy  +  ...  Gp  Sp  +  G)=  Gy*  VCGi)  +  . . .  +  Gp  *  V(Sp)  +  G  (73) 

We  will  call  such  af,a  value  function.  Notice  that  for  any  value  function  V  and  any  interval  G, 
V(G)  —  G. 

Lemma  1.  If  V  associates  with  each  symbol  a  singleton  interval  17(5)  =  [s,  s]  then  for  any 
[A,  A']  E  A,  and /  E  0. 


cl  E  V (A),  cl'  E  V (A)  =$  cl  +  cl'  G  V (A  +  A?) 

(74) 

a  E  V (A),  a'  EV (A')  =>  a  —  a'  E  17 (A  — 

A') 

(75) 

aE  V(A),iE0  =>  ai  E  V(A*I) 

(76) 

a  E  17(A),  i  EJ  ^  a/i  E  V(A/I) 

(77) 

a  E  V(A),Si  E  V(A ')  =>  a  E  V(Subst[A,Si 

-»A') 

(78) 

Proof:  Let 

A  —  Gy  Sy  +  ...  Gp  Sp  +  G 

(79) 

A'  =  G'ySy  +  ...  G'pSp  +  G' 

(80) 

a  =  gySy  +  ...  +  gpSp  +  g  with  gt  E  Gt, g  E  G 

(81) 

a'  =  g'ySy  +  . . .  +  g'psp  +  g'  with  g\  E  Gu  g'  E  G 

(82) 
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(83) 


Then,  using  distributivity,  we  get 

a  +  a'  =  (flfjl  +  g'Js!  +  ...  +  (jgp  +  g'p)sp  +  (g  +  g') 

£  (Gi  +  +  ...  +  ( Gp  +  G' p)Sp  +  (G  +  G*) 

£  V (A  +  A) 

The  other  cases  are  similar.  Notice  that  the  first  two  implications  do  not  necessarily  hold  if 
V  (S’)  is  not  a  singleton  interval  as  interval  multiplication  is  not  distributive  over  interval  addition. 
QED 

4.1.7. 2.  Error  Set  Abstractions 


For  each  variable  x,  we  will  define  F(x)  £  c/Z  and  then  show  how  it  represents  an  abstraction  of 
the  set  of  errors  for  variable  x.  More  precisely,  F(x)  will  represent  a  symbolic  combination  of 
the  inputs  and  outputs. 

Definition  of  E(x).  For  constants  c,  let  F(c)  be  a  constant  interval  that  includes  all  possible 
error  values  for  c. 

For  any  input  ut,  let  t/j  £  S  and  let  E(ui)  —  Ui  —  [1,1]  1/j  £  c/Z. 

Similarly,  for  output  y£,  let  Yt  £  S  and  E(yi)  —  Yt  —  [1,1] Y[  £  c/Z. 

Suppose  that  for  all  x  that  were  assigned  a  value  before  v  in  the  filter  code,  E  (x)  £  c/Z  was 
defined.  Then  the  definition  for  E(y')  depends  on  the  statement  that  assigns  it  a  value. 

v  —  w  ■ 

E(y)  =  E(w )  (84) 

v  —  c- 


E{y )  =  E(c) 

V  =  W  +  Z  : 

E(v)  =  R(v)  *  [-em,em]  +  (E(w)  +  E(z))  *  [l-em,l  +  £m] 
v  —  w  —  z  ■ 

E(v)  =  R(v)  *  [-em,em\  +  (E(w)  -  F(z))  *  [1  -  em,  1  +  em  ] 

V  —  c  *  w  ■ 

E(v)  =  R(v)  *  [-em,em]  +  (f(c)  *  R(w)  +  (R(c)  +  E(c))  *  F(w))  *  [l-em,l  +  em] 

V  —  w  /  c  : 

„ ,  .  ~R(v)  E(c )  +  R(w)  *  [~em,  em]  +  E(w)*[  1  -  em,  l  +  em] 

^  = - m  +  m - 


(85) 

(86) 

(87) 

(88) 

(89) 
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Lemma  2.  Let  Cj(it)  and  eL (y)  be  the  error  for  input  i  and  output  i  respectively  and  V  the  value 
function  with  F(l/j)  =  [fj(it),  Cj(it)]  and  V(Yi)  —  [^(y),  £j(y )  ].  Then  for  any  variable  x 


Vn  >0.fi(i)£f(f(x)).  (90) 

Proof:  From  the  definition  of  E  (x),  the  statement  is  trivially  true  for  the  state  variables  and 
constants.  Assume  that  the  statement  is  true  for  the  variables  given  a  value  before  v.  Then  the 
first  two  cases  are  obvious. 

Case  v  =  w  +  z: 

fn(y^  —  Yn(y) 

=  ( fn(w )  +  /n(z))  (1  +  5)  -  (rn(w)  +  rn(z))  where  \S\  <  em 
=  0n(w )  +  e„(w))+  (rn(z)  +  en(z)))  (1  +  5)  -  (rn(w)  +  rn(z))  (91) 

=  (r„(w)  +  rn(z))  8  +  (en(w)  +  e„(z))  (1  +  5) 

=  rn(i?)  8  +  (en(w)  +  en(z))  (1  +  5) 

then, 

rn(u)  G  R(v)  =  V(R(v)), 

5  £  [  £jn<  ^ml  —  ^([— ^m>  ^mJ)> 

1  +  5  G  [1  -  em,  1  +  em]  =  F([l  -  e™,  1  +  Cn]  ) 

and  from  the  inductive  hypothesis  we  have,  en(w)  G  V (E(w))  and  en (z)  G  F (E(z)).  Therefore, 
by  Lemma  1 , 

/n(v)  -  e  *  [-em-  fm]  +  (  E (w)  +  £(z)  )  *  [1  —  em,  1  +  ej) 

= 

Case  v  =  w  —  z:  similar  to  previous  case 

Case  u  =  c  *  w: 

/nW-  rn(l?) 

=  /(c)  /n(w)  (1  +  8)  -  r(c)  rn(w)  where  |  5|  < 

=  (r(c)  +  e(c))  (rn(w)  +  e„(w))  (1  +  5)  -  r(c)rn(w) 

=  r(c)  rn(w)  5  +  (e(c)rn(w)  +  r(c)  en(w)  +  e(c)  en(w))  (1  +  5) 

=  rn(v)  5  +  (e(c)rn(w)  +  (r(c)  +  e(c))  en(w))  (1  +  5) 

G  F(/?(v)  *  [-em,em]  +  (E(c)  *  R(w)  +  (R(c)  +  E(cj)  *  E(w ))  *  [1  -  em,l  +  em]  ) 

-  F(E(iO) 

Case  v  —  w  /  c\ 
fn(v)~  fit 
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=  (/nO)  //(c))  (1  +  5)  -  rn(w)  /r(c)  where  |  5|  <  em 
(rn(w)  +  en(w))  (1  +  5)  _  rn(w) 
r(c)  +  e(c)  r(c) 

=  (rn(w)  +  en(w))  (1  +  5)  -  rn(w)  (r(c)  +  e(c))  /r(c) 

r(c)  +  e(c) 

(rn(w)  r(c)  -  rn(w)  r(c)  -  rn(w)  e(c))  /r(c)  +  rn(w)  5  +  en(w)  (1  +  5) 

r(c)  +  e(c) 

-  rn(w )  e(c)  /  r(c)  +  rn(w)  5  +  en(w)  (1  +  5) 


(93) 


r(c)  +  e(c) 

-  flQO  *  g(c)  +  fl(w)  *  [-e^grn]  +  E(w)  *  [1  -  1  +  em\ 

V  V  n  f  i  T7  f  ) 


R(c)  +  £(c) 


=  E(E(v)). 

QED 


Notice  that,  for  any  x,  E(x)  is  either  a  symbol  (in  the  case  when  x  is  an  input  or  a  past  output), 
a  constant  interval  (in  the  case  x  is  a  constant  value)  or  is  defined  as  a  linear  expression  over 


other  error  set  abstractions.  Therefore  any  E(x)  can  be  expressed  as  a  linear  combination  of 
input  symbols  Lf  and  past  output  symbols  V). 

Suppose  that  for  the  output  variable  y,  E  (y)  G  c/Z  is 

£(y)  =  ^n+l+  ■  ■  ■  +  -dw  Un-N+ 1  +  Tn  +  . . .  +  Y-n-N+l  +  P  (94) 

Without  loss  of  generality  we  can  assume  P  is  symmetric  about  0.  Indeed,  we  can  always  enlarge 
P  into  a  symmetric  interval  and  the  effect  will  be  a  larger  E  (y) .  We  can  treat  P  similar  to  an 
input  error  at  iteration  n  +  1  by  introducing  a  symbolic  variable  P-n+i  ■ 

E(y)  =  A0Un+1+ ...  + AN  Un-N+1  +  B1Yn  +  ...  +  BNYn_N+1  + [1,1]  Pn+1  (95) 
We  will  define  recursively  En(y )  G  c/Z  to  be  the  unfolded  version  of  E(y)  of  “depth”  n. 

Ek(y)  =  [0,0]  if  k  =  0 . N-  1  (96) 

£fc+i(y)  =  Subst(E(y),  Yn  ->  Ek(y),...,Yn_N+1  ->  Efc_w+1(y))  for  k  >  N  (97) 
With  these  definitions, 

£n+l(y)  =  C0  t/n+1+...+ Ctf  t/i  +  ^0  ^n+l  +  ■  ■  ■  +  DnPl  (98) 

where  Cn,  Dn are  intervals  and  for  n  >  iV, 

Cn+i  =  5iCn  +  ...  +  BNCn_N+1  (99) 

Pn+l  =  P\Pn  +  ■■■+  PnPu-N+1  (10°) 
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4. 1.7. 3.  Error  Bounds  for  First-Order  Linear  Filters 

For  a  first  order  filter,  the  recurrence  equations  are 


Co  =  A0  (101) 

Ci  =  A1  +  B±A0  (102) 

Cn+i  =  B1Cn  (103) 

Co  =  [1,1]  (104) 

T>i  =  fi^U]  =  B1  (105) 

Dn+ 1  =  BxDn  (106) 

which  have  the  solution 

Cn+1  =  C.Bf  (107) 

Dn+1  =  DxBl  =  (1°8) 


Let  V  the  value  function  for  which  V  (Lf )  =  EV (iq)  (the  set  of  error  values  for  iq)  and  V (fij)  = 
P.  Then 


V(En+1(y))  =  C0  EV (un+1)  +  C1EV(un)  +  ...C1B^EV(u1)  + 
D0  EV(pn+1 )  +  D1  EV(pn )  +  . . . 

For  an  interval  /  =  [a,  b],  let  max  |/|  =  max{  |x|  |  x  G  [a,  b]}  —  max(\a\, \b\). 
As  EV  (it)  and  P  are  symmetric  about  0, 


max  \V(En+1(y))\  <  max  |C0|  max(EV(u))  +  max(P)  + 

max  |n  max(EV(u ))  +  . . .  +  max  \CX \  (max  |fii|)n  max(EV (it))  + 
max  \D1\  max(P)(  1  +  . . .  +  (max  |fiil)n) 

=  max  \C0\  max(EV (u))  +  max(fi)  + 

(max  m  max(fiF(it))  +  max  |fix|  max(fi))(l  +  . . .  +  (max  IfljJ)”) 

For  example,  let’s  reconsider  the  lead  lag  filter  described  by 

y0  =  0  (110) 

9  7  11 

Vn+l  ~  TtT^n+l  i  0  tin  +  1  0  In  (111) 


with  u0  —  0  and  for  any  t  >  0,  iq  G  [—1, 1], 

The  analysis  produces  the  following  recurrence  equation  for  the  accumulated  error  after  n  +  1 
iterations 
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En+iiy)  =  A0EV(un+ 1)  +  A1EV(un)  +  5i£F(yn)  +  EV(pn+1 ) 


(112) 


where 


Therefore, 


rO. 69230720700609994212820001624950,1 
0  ~~  I-0.69230816632349415561292123758108J 

A1  =  [-0.53846195604689442353836541364470, 

-0.53846117730685048103835362264472] 

B1  =  [0.84615337425443068269699198003444, 

0.84615438577002839040341485924719] 

EVipn)  =  P  =  [-0.00000073359592004213967781243823, 

0.00000073359592004213967781243823] 

\V(En+1(y))\  <  Aerror  3”  M error  (1  +  ■  ■  ■  +  ) 


where 


C0  =  A0  =  [0.69230720700609994212820001624950, 
0.69230816632349415561292123758108] 


(113) 

(114) 

(115) 

(116) 

(117) 

(118) 


C1=A1  +  B1C0  =  [0.04733612318197767638922552809080, 

0.04733841393218037023261187048824]  (119) 


Aerror  =  max\C0\max(EV  (u))  +  max(P)  =  0.00000081612548469976751733157194 
terror  =  max|Ci \max(EV(u))  +  7nax|51|max(P)  =  0.00000062637858381997158563381300 
p  =  maxIBil  =  0.84615438577002839040341485924719 

4.1. 7.4.  Error  Bounds  for  Second-Order  Linear  Filters 

For  a  second-order  filter,  the  recurrence  equations  are 


II 

O 

(120) 

Cl 

—  A1  +  B1C0 

(121) 

C2  ~  A2  +  B1C1  +  B2C0 

(122) 

C-n+i  —  Bi  Cn  +  B2Cn- 1 

(123) 

D0  =  [1,1] 

(124) 

Di 

=  fii[l,l]  =  Si 

(125) 

D2  =  B1D1  +  S2[l,l]  =  Bl  +  B2 

(126) 

Dn+ 1  —  D1Dn  +  B2Dn_  1 

(127) 

Approved  for  Public  Release;  Distribution  Unlimited. 
38 


(128) 


Consider  all  the  solutions  rt  <  r2  of  the  equation 

x2  -  bxx  +  b2 

with  bx  E  B±  and  b2  E  B2.  Then 

Bi  —  -J  B?  +  4fl2 
«*=  -  2 

Bi  +  - J  B?  +  4B2 

R2  =  — — - - 

2  2 

and  therefore 

{rx  |  rx  the  smaller  solution  ofx2  =  bxx  +  b2  for  bt  E  Bt}  Q 
{r2  |  r2  the  larger  solution  ofx2  =  brx  +  b2  for  bt  E  SJ  Q  R2 


(129) 

(130) 


In  the  case  when  Bjand  R2  do  not  intersect,  we  can  define  a  few  constants  that  are  needed  in 
what  follows 


maxtRi)  max(R2  —  R±)  +  max(R2 )  max(R2  —  R2) 
^  ~  min(R2  -  Rx) 

max)!?!)  max(R1  —  Bx)  +  max(R2 )  max(R2  —  R^ 
~~  min(R2  -  /?,) 

/?  =  max  (/?!,/?2) 


M 

N 


max  1 0^2  —  C2|  max  |C2  —  C1R1\ 

max  ( - - - - — , - - - — ) 

min{R2  —  Bx)  min{R2  —  Rx) 

max  \D1R2  —  D2\  max  \D2  —  D1R1\ 

max( - - - - — , - - - - — ) 

min(R2  —  Rx)  min{R2  —  Rx) 


(131) 

(132) 

(133) 

(134) 

(135) 


Theorem:  If  0  <  min^R^  and  max(R1)  <  min(R2 )  then  for  any  n  >  1, 

\En+2(y)\  <  [max|C0|  +  max|C!|  +  max|C2|  +  2Mmax(R2)  (1  +  . . .  +  fi11-1)]  max(EV(u))  + 
[max  |D0|  +  max  \D1\  +  max|D2|+  2  A/ rnax(i?2)(l  +  ■  ■  ■  + /?n_1)]  wax)?)  (135) 

The  theorem  implies  that 

max  \  V(En+1(y))\  <  Aerror  +  Merror  (1  +  •••  +  /?”)  (137) 


where 

Aerror  —  [max|C0|  +  max | Cx\  +  max|C2|]  max(£V(it))  + 

[max  \D0  \  +  max  |DX|  +  max  \D2\]  max(P) 

Merror  =  2  Mmax(R2 )  max(EV  (u))  +  2  N  max(R2 )  max(P) 


(138) 

(139) 
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Notations:  Let  c0,  c1, . ..  be  such  that  c0  G  C0,  c±  G  C±  and  V  n  >  0.  cn+1  =  bn+11  cn  + 
bn+12  cn-i  with  bn+l  i  G  flj.  For  each  n  >  0,  let  rn+11<  rn+12  be  the  solutions  of 

x2  =  bn+lilx  +  bn+ii2  (140) 


and  let  sn+  l  and  tn+1  be  the  solutions  of 

cn- 1  =  sn+l  "b  tn+1 
cn  ~  sn+l  rn+l,l  +  tn+ 1  rn+ 1,2 


(141) 

(142) 


Lemma  3:  If  0  <  min^R^  and  max^R^  <  min(R2 )  then 
max  (\s3\,\t3\)  <  M 

max  (|5n+2|,  |tn+2|)  <  max  (|sn+1|,  |tn+1|)  *  /? 

Proof:  s3 

and  t3  are  the  solutions  of  the  system 

Cl  =  s3  +  t3 

(143) 

c2  =  s3  r31  +  t3r3]2 

(144) 

Therefore 

Cir3,2  ~  c2 

S3  = 

r3,2  r3,l 

(145) 

c2  ~  Clr3,l 

to  — 

r3,2  —  r3,l 

(146) 

which  imply  the  first  inequality  of  Lemma  3.  Now,  for  any  n  >  1, 


cn+ 1  _  ^n+1,1  cn  +  ^n+1,2  cn-l  —  ^n+1,1  (sn+l  rn+l,l  +  tn+1  rn+l,2)  +  ^n+1,2  (sn+l 
—  sn+ 1  (^n+1,1  fii+1,1  +  ^n+1,2)  +  tn+i(hn+i  1  rn+12  +  bn+12) 

—  Ni+1  fii+1,1  +  tn+1  rn+12 
From  the  definitions  of  rn+l  i  and  rn+2i, 

cn  ~  5n+l  fii+1,1  +  tn+1  rn+l,2  =  sn+ 2  +  tn+2 
cn+l  ~  sn+ 1  rn+l,l  +  tn+1  T^+1  2  =  Sn+2  Tn+21  +  tn+2  rn+2  x 


+  tn+i) 

(147) 


(148) 

(149) 
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and  by  solving  this  system  for  sn+2  and  tn+2  we  get 


Vn  >  0.sn+2 
Vn  >  0.tn+2 


—  sn+ 1 

—  sn+l 


rn+l,l\rn+2,2  rn+ 1,1, 
rn+ 2,2  ~  rn+ 2,1 
rn+l,l(  rn+l,l  —  rn+ 2,1 
rn+2,2  ~  rn+ 2,1 


+  £ji+1 
+  Li+1 


rn+l,2(  rn+2,2  ~  rn+l,2) 
rn+2,2  ~  rn+ 2,1 
rn+l,2(  rn+l,2  ~  rn+2,l) 
rn+ 2,2  ~  rn+ 2,1 


(150) 

(151) 


therefore 


l^n+2 I  —  lsn+ll 
I  Li+2 I  —  l^n+2 I 


max(/?1)  max{R2  —  Rx) 

+  iLi+il 

max{R2 )  max(R2  —  R2 ) 

(152) 

min(R2  —  R±) 

min{R2  —  Rx) 

max^Ri)  max(R1  —  Rx) 

+  iLi+ll 

max(R2 )  max{R2  —  R x) 

(153) 

min(R2  —  Rx) 

min{R2  —  Rx) 

and  thus 


|sn+2|  <max  (|sn+1|,|tn+1|)  /?!  (154) 

|tn+2|  <  max  (|sn+1|,  |tn+1|)  /?2  (155) 

which  imply  the  second  inequality  of  Lemma  3.  QED 

Proof  of  Theorem:  From  Lemma  3, 

max(\sn+2\,\tn+2\)  <  max  (\sn+1\,\tn+1\)P  <  M  /?n_1  (156) 


Thus 


which  implies 


I  cn+2 1  <M/3n  1  rx  +  M/?n  1  r2  <  2 Mmax(R2)  /3n  1 
max|Cn+2|  <2 Mmax(l?2)  /?n'1 


Similarly,  for  the  sequence  Dn, 

max  \Dn+2\  <  2Nmax(R2)  /?n_1 


(157) 

(158) 


(159) 


Thus 

£n+2(y)  e  C0^(un+1) +  ...+  CnEV(u x)  +  D0  EV(pn+1)  +  . . .  +  DnEVip-i) 
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\En+2(y)\  <  max  |  C0|  max(EV (it))  +  ...  +  max  \  Cn  \  max(EV(u ))  + 

max  |  D0\max(P)+  . . .  +  max  \  Dn\max{P )  (160) 

=  ( max\C0\  +  maxlCj  +  max\C2\  +  2Mmax(/?2)(l  +  — 1-  /?n_1)  max(EV(u)')  + 

max\D0\  +max|D1|  +  max\D2\  +  2Nmax(R2)  (1  +  — 1 - /?n_1)  max(P) 


QED 

Example:  consider  again  the  second-order  filter  (from  Section  4.2.4) 

yo  =  7i  =  0  (161) 

(162) 

592841  i  1185682  _  592841  _  126814318  51175081 

yn+1  ~  78010601  Un+1  +  78010601  Un  +  78010601  Un~x  +  78010601  yn  ~  78010601  yn_1 

where  u0  =  ux  =  0  and  for  any  i  >  0,  iq  £  [—100, 100]. 

The  analysis  produces  the  following  recurrence  equation  for  the  accumulated  error  of  the  output 
£n+i(y)  =  A0EV(Un+1)  +A1EV(un)  +A2EV(un_1)  +  B-^EViyn)  +  B2EV(yn_1 )  +  EV(pn+1 ) 

(163) 


where 

A0  =  [0.00759948777424905218315802011899, 
0.00759950078268356597649013476543] 

Ax  =  [0.01519897917221431452601296749856, 
0.01519900518909170492103914110864] 
d2  =  [0.00759948868017869869646495929948, 
0.00759949987675258472756483898753] 

S,  =  [1.62560296051380471713218194062187, 
1.62560458404225805892680063197886] 

B2  =  [-0.65600182747582427514452884015165, 

-0.65600 12873574070 1 548264934562427] 

P  = [-0.000039589480719188901 18042099520, 
0.0000395894807 191889011 8042099520] 

C0=A0  =  [0.00759948777424905218315802011899, 
0.00759950078268356597649013476543] 
C1=A1+  BtC  o  =  [0.02755272899642203819834661633796, 
0.02755278849785483774577685006890] 

c2  =a2  +b1c1  +  b2c0 

=  [0.04740400010565253933432332363906, 
0.047404165398847956558342881 171 12] 

D0  =  [hi] 
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D  i  =  B1  [1,1]  =  B1  =  [1.62560296051380471713218194062187, 
1.62560458404225805892680063197886] 
D2  =  B1D1  +  B2  [1,1]  =  B\  +  B2 

=  [1.98658315775542226318354179506256, 
1.98658897630179582912383485351786] 

D  _  Bt  -  y]Bl  +4  B2 
ri-  2 

=  [0.744637864 1 747600052 1570285455046, 

0.7446523 1 893245232206672049472825] 

B1  +  Vfii2  +4  B2 


=  [0.88095145334557906596277079157211, 
0.88096590810327138281378843174991] 

max(l?1)  max(R2  —  Rx)  +  max(R2 )  max(R2  —  R2 ) 
min(R2  —  Rx) 

=  0.74490369020164437255103955622501 


P 2  = 


max(/?1)  max(R1  —  Rx)  +  max(R2 )  max(R2  —  Rx) 
min(R2  —  Rx) 


=  0.88123173566310767125443267351180 
0=  maxiP^PJ  =  0.88122124050850451432751993407794 


max\ClR2  —  C 2|  max \C2  —  C1R1\ 

M  =  TTICLX  ( -  - 

min(R2  —  R 1)  '  min(R2  —  R 1) 

=  0.19726728451018537213720993866727 
max  \D1  R2  —  D2\}{ }  max  \D2 


) 


N  —  max(- 


D 1  R  1| 


min{R2  —  R 1) 


min(R  2  —  R 1) 


■) 


=  5.69411877140933861544292553129993 
terror  =  0.00018357849766142422504207973351 
terror  =  0.00040133074049677426939408818760 


(164) 

(165) 

(166) 

(167) 

(168) 

(169) 

(170) 

(171) 

(172) 


4.1.8.  Accumulating  Integrators:  Error  Range  Functions 

In  contrast  to  digital  filters,  accumulating  integrators  have  no  discounting  mechanism  for  past 
inputs.  Instead,  the  output  range  intentionally  grows  without  bound.  An  analysis  of  the 
floating-point  errors  can  only  give  bounds  as  a  function  of  the  number  of  iterations.  The  analyst 
can  then  decide  how  many  iterations  to  allow  before  the  error  becomes  unacceptable. 

Let 
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(173) 

(174) 


*o  =  0 

*n+ 1  —  *n  4”  c- 


be  the  recursive  equation  of  an  integrator.  In  this  case,  we  are  interested  in  the  range  of  values 
after  at  most  n  iterations,  which  is  [0,  nc]  or  [ nc ,  0]  depending  whether  c  is  a  positive  or  negative 
constant. 

To  find  the  error  range  after  n  iterations  or  less,  we  cannot  use  the  same  expressions  as  for  the 
analysis  of  filters.  In  that  case,  we  used  the  fact  that  all  the  variables  involved  have  finite  ranges 
which  we  computed  in  advance.  However,  we  can  adapt  the  method  by  changing  the  symbolic 
interval  combinations  used  for  filters  with  symbolic  combinations  over  linear  expressions  in  n, 
the  number  of  iterations. 

The  real  value  after  n  iterations  is  Rn(x)  =  nc.  In  the  simple  case  when  c  is  found  to  be 
representable  then  the  error  found  is  0.  Otherwise,  the  analysis  finds  the  following  recurrence 
relation. 


En+ iM  =  En(x)B1  +  Qn  +  P 

(175) 

Bi=  [1  -  em,  1  +  em] 

(176) 

Q  —  C  *  [  €m,  Cjjj] 

(177) 

P  =  E(c)  *  [1  -  em,  1  +  em] 

(178) 

from  which  we  get  the  closed  fonn 

\V(En+1(x))\  <  max  \Q\  (n  +  (n  —  1)/?  +  •••  +  /?n)  +  max  |P|  (1  +  b (3  +  ■■■  +  /?n) 


where  /?  =  1  +  em.  Notice  that 
n  +  (n  —  1)/?  +  •••  +  /?n 

=  (1  +  b/3  +  -  +  /?n)  +  (1  +  bp  +  -  +  Z?”-1)  +  ...  +  1 


(179) 


on+l  _  on  _  o  _  i 

P  +  ^ - —  +  •"  +P 


/?-  1  /?-  1 
/?  +  -  +  /?n+1 


P-  1 


n 


P-1 


B  n 

(1  +-  +  /?n)  - 


P-1 


/?-! 


(180) 
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and  therefore 


|K(£n+1(x))|  <  max\Q\-^-j(l  +/?  +  •••  +  /?n)  -  max\Q\  + 
max  |P|  (1  +  /?  +  — I-  /?n) 

or 

max  \V (En+1(y))\  < 

Aerror  +  Merror(l  +  (3  +  •••  +  (3  )  +  Nerrorn 

where 

a  —  n 

nerror  u 

/? 

terror  =  max|Q|  — -  +  max  \P\  =  (|c|  +  E(c))(l  +  6m) 

max  1 0 1 
N  - - — 

1V error  ^  ^ 

Example:  Consider  the  integrator  described  by 

(181) 

(182) 


*o  =  0 

xn+i  xn  +  0.1 


with  the  output  variable  given  by  yn  =  xn  *  3.2.  This  is  the  essence  of  the  integrator  that 
caused  the  Patriot  Missile  failure.  The  error  bound  for  the  integrator  is  then 

max  |P(£n+1(x))|  < 

Ax, error  ~b  ^x, error  (1  H —  +  /?n)  +  Nx  error  n 

where 


A 


x,  error 


M. 


x, error 


N. 


x, error 


=  o 

=  (|c|  +  £(c))(  l  +  em)  =  0.10000001788139414315992326010019 
=  -|c|  =  -0.1 


The  error  bound  on  the  output  is 

En(y)  =  Rn(y)  *  +  [E (3.2)  *  £n(x)  +  (3.2  +  £(3.2))  *  £n(x)]  *[l-em,l  +  em] 

(183) 


and  thus 


max  |P(£n+i(y))l  ^  0.32  emn 

'  +  [£( 3.2)  *0.1  n  +  (3.2  +  £(3.2))  *  (M*,error(  1  +•••  +  /?")-  0.1  n)]  (1  +  em) 
Ay, error  T  My  error  (1  T  T  (3  )  +  Ny, error  (184) 
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where 


a  —  n 

riy  error  u 

My  error  ~  (3.2  +  E{  3.2))(1  +  em)Mx 

error 

=  0.32000011444093274803971696172939  (185) 

Ny, error  =  0-32  6m  +  0.1  E( 3.2)(1  +  em)  -  (3.2  +  £'(3.2))0.1  (1  +  em) 

=  0.32  em  -  0.32  (1  +  em)  =  -  0.32 


4.2.  Summary  of  Filter  Test  Suite  Results 

The  development  of  the  first  year  analysis  of  Linear  Digital  Filters  was  guided  by,  tuned,  and 
ultimately  tested  for  round  trip  integration  by  a  suite  of  test  cases  designed  by  the  Honeywell 
team.  These  consisted  of  Simulink  models  of  representative  filter  applications,  their  generated  C 
implementations  and  their  generated  JSON  files  describing  model-level  properties  of  the 
generated  code  to  be  exploited  by  CodeHawk.  This  section  summarizes  the  results  of  each  of 
these  test  cases. 

General  Significance  and  Organization  of  the  Results. 

As  discussed  in  Section  3.1,  the  goal  of  the  combined  model-level  and  code-level  analysis  is  to 
derive  precise,  consistent  bounds  for  the  range  and  numerical  error  of  the  filter  outputs.  Both  the 
range  and  error  need  to  be  bounded  in  order  to  guarantee  stability  of  the  control  algorithm.  As  we 
have  pointed  out  earlier  in  this  report,  most  static  analysis  tools  are  not  able  to  compute  such  a 
bound  and  can  thus  report  high  bounds.  This  forces  current  practices  to  rely  upon  informal, 
empirical  manual  analysis/reviews  and  simulations  of  the  code.  Our  results  demonstrate  that  our 
combined  analysis  technique  proves  precise,  conservative  bounds  for  both  variable  range  and 
numerical  error  -  representing  a  definite  advancement  in  the  state-of-art  of  static  analysis 
approaches  for  numerical  algorithms. 

Range  Bounds  Results:  The  range  at  the  output  of  the  filter  is  theoretically  expected  to  be 
bounded.  The  bounds  should  be  stable  over  time  -  often  proportionally  related  to  the  input  range 
bound  based  upon  the  filter  coefficients  which  are  derived  from  the  time  constants  (Tn  and  Td)  of 
the  filter  and  the  sample  period  (Ts).  As  presented  in  the  following  subsections,  the  results 
returned  by  the  CodeHawk  analysis  of  filter  range  bounds  on  he  source  code  confirm  the 
theoretically  expected  properties.  Furthermore  the  bounds  are  precise  (tight),  yet  conservative  - 
i.e.,  representing  the  worst  case  if  value  at  the  filter  input  is  varied  within  its  range  arbitrarily 
over  time. 

Numerical  Error  Bounds  Results:  For  a  properly  designed  filter,  the  numerical  error  at  the 
output  of  the  filter  is  expected  to  be  bounded  to  a  small  value  (less  than  1%  of  the  range)  in  order 
to  achieve  control  stability.  The  error  bounds  should  be  convergent  over  time.  The  ideal  dynamic 
scenario  is  that  the  numerical  error  should  approach  0  as  time  approaches  infinity  if  the  filter 
input  is  held  at  a  constant  value.  This  ideal  dynamic  scenario,  however,  is  hard  to  achieve  in  a 
static  analysis  framework.  Thus,  industry  practice  relies  on  empirically  testing  the  tolerance  for 
error  to  be  within  1%  of  the  variable  value.  Our  analysis  proves  the  error  bound  to  be  well  within 
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this  tolerance  even  for  2nd  order  (quadratic)  filter;  furthermore  the  error  bound  remains  constant 
with  time. 

4.2.1.  Filter  Example  using  a  Resettable  Lead  Lag  Filter 


Figure  13  Resettable  Lead  Lag  Filter  -  Test  Model  Diagram 

Figure  13  shows  ‘errorcomp’,  a  model  with  two  resettable  lead  lag  filters  in  series,  creating  a  2nd 
order  filter.  The  next  subsection  describes  the  transfer  function  and  difference  equation 
derivation  for  this  class  of  filters.  Following  subsection  shows  the  analysis  results  assuming  the 
filter  runs  for  72,000  seconds  at  40  Hz  (25  ms  period). 

4.2. 1.1.  Resettable  Lead  Lag  Filter  Transfer  Function 

Figure  14  shows  the  Resettable  Lead  Lag  filter  as  it  appears  as  a  block  in  MATLAB  Simulink. 
This  class  of  filters  has  two  time  constants,  Tn  and  Td,  that  are  both  constants  for  a  particular 
filter  instance.  There  are  additional  reset  (IC)  and  reset  value  (ICV)  inputs  that  allow  the  filter  to 
be  reset. 


"t 


Figure  14  Resettable  Lead  Lag  Filter 


In 


ICV 

IC 


Tn  s  +1 


Td  s  +1 


istLeadLag 
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The  following  is  the  transfer  function  derivation  for  this  filter  class: 


Reset  Mode  (when  IC  —  1):  y  —  ICV 
Filter  Mode  Transfer  Function  (when  IC  —  0): 


General  form: 


T  s+1 

Continuous  domain:  H(s )  =  — — 

V  '  Tds+ 1 

Coefficient  definitions  (all  coefficients  are  constant) 

2  Td 


Discrete  domain:  H(z )  = 


/C3+/C4Z 

k\+k-2Z~ 


2  Td 

fei  =  1  +  — - 


k7  =  1  - 


ko  —  1  + 


2r„ 


k,  =  1  - 


2T„ 


Difference  equation: 


yn  =  (fc3“n  +  k4“n-l  -  *2yn-l)/*l 
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(186) 


4.2.I.2.  Analysis  results  after  72,000  Seconds  (20  hrs) 


System  OperTime  (secs):  72000  Model  Period  (ms):  25  Total  Periods:  2880000 

#  Symbol 

Name 

Normal 

Range  Min 

Normal 

Range  Max 

Range 

Constraint 
Min  (if  any) 

Range 

Constraint 
Max  (if  any) 

Data  Type 

Worst  Case 

Numerical  Error 

Bounds 

CMP_ERROR 

-11.65675964 

11.65675964 

-11.65675964 

11.65675964 

Langlndep_Float32 

eB[-+1.33e-003] 

#  Filters 

cmp3 

-3.950617284 

3.950617284 

-3.950617284 

3.950617284 

Langlndep_Float32 

eB[-+2.29e-004] 

cmp4 

-11.65675964 

11.65675964 

-11.65675964 

11.65675964 

Langlndep_Float32 

eB[-+1.33e-003] 

Figure  15  Resettable  Lead  Lag  Filter  -  72000  Seconds 

In  Figure  15,  variable  CMPERROR  shows  the  analysis  for  the  final  output  of  the  model  in 
Figure  13  for  a  run  time  of  72000  seconds.  Cmp3  and  cmp4  show  the  analysis  for  the  outputs  of 
filters  cmp3  and  cmp4  respectively.  The  output  of  cmp4  matches  the  values  of  CMP  ERROR 
because  there  is  no  intervening  functionality.  For  the  filter  cmp3,  the  normal  (operational)  range 
of  the  filter  input  is  [-2,  2]  which  also  happens  to  be  a  hard  range  constraint  since  there  is  a  limit 
block  before  cmp3,  constraining  the  range  not  to  exceed  [-2,  2]  even  if  abnormal  values  are 
present  at  the  input  variable  of  the  model  named  ERROR.  Note  that  the  range  bounds  reported  in 
Figure  15  are  convergent  over  time;  i.e.  filter  output  range  is  bounded  by  the  values  shown  when 
input  value  varies  arbitrarily  within  its  range  over  infinite  time.  The  range  bounds  of  [-3.9506, 
3.9506]  for  cmp3  output  are  precise  and  compliant  with  the  theoretically  expected  bounds  for  the 
specific  values  of  filter  coefficients  for  cmp3.  Likewise  for  cmp4. 

The  worst  case  numerical  error  bound  (WCEB)  for  cmp3  output,  after  72000  seconds  of 
operation,  is  [-2.29e-004,  -2.29e-004]  over  the  range  of  [-3.9506,  3.9506]  -  this  is  two  orders  of 
magnitude  less  than  the  empirical  practice  tolerance  of  1%. 
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4.2.2.  Filter  Example  using  a  Variable  Lag  Filter 


Figure  16  Variable  Lag  Filter  -  Test  Model  Diagram 

Figure  16  shows  ‘highlimsetpoint’,  a  model  containing  a  variable  lag  filter  with  a  fixed  tau. 
The  next  subsection  describes  the  transfer  function  and  difference  equation  derivation  for  this 
class  of  filters.  Following  subsection  shows  the  analysis  results  assuming  the  filter  runs  for 
72,000  seconds  at  10  Hz  (100  ms  period). 

4.2.2. 1.  Variable  Lag  Filter  Transfer  Function 

Figure  17  shows  the  Variable  Lag  filter  as  it  appears  as  a  block  in  MATLAB  Simulink.  This 
class  of  filters  has  a  variable  tau,  x,  that  determines  the  lag  according  to  the  transfer  function. 
There  are  additional  reset  (IC)  and  reset  value  (ICV)  inputs  that  allow  the  filter  to  be  reset. 

> 

> 

> 

> 


In 

ICV 

lc  T5 

7 

+-  1 

varLag 


Figure  17  Variable  Lag  Filter 
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The  following  is  the  transfer  function  derivation  for  this  filter  class: 

Reset  Mode  (when  IC  =  1):  y  =  ICV 

Filter  Mode  Transfer  Function  (when  IC  =  0): 

General  form: 

1  1 +z~1 

Continuous  domain:  H(s )  = -  Discrete  domain:  H(z)  = - T 

v  '  ts+1  v  J  k1+k2z~1 

Coefficient  definitions  (r  is  given  as  an  input) 


Difference  equation: 

yn  =  (un  +  un_!  -  fezyn-O/fei  ( 1 87) 

4. 2.2.2.  Analysis  results  after  72,000  Seconds  (20  hrs) 


System  OperTime  (secs):  72000  Model  Period  (ms):  lOO  Total  Periods:  720000 

#  Symbol  Name 

Normal 

Range  Min 

Normal 

Range  Max 

Range 

Constrai 

nt  Min 

(if  any) 

Range 

Constrai 

nt  Max 

(if  any) 

Data  Type 

Worst  Case 

Numerical  Error 

Bounds 

CLOSED_LOOP 

O 

1 

Langlndep_Bool 

e  B  [ -+0 .  OO  e +000 ] 

HI  LIMC 

1.239987848 

130.6000137 

_ o 

132 

Langlndep_Float32 

eB[-+4.53e-004] 

#  Filters 

varLag 

O 

0.980000196 

o 

1 

Langlndep_Float32 

eB[-+2.34e-007] 

Figure  18  Variable  Lag  Filter  -  72000  Seconds 

In  the  above  figure,  CLOSED  LOOP  is  Boolean,  and  thus  incurs  no  numerical  error.  HILIMC 
shows  the  analysis  for  the  final  output  of  the  model  in  Figure  16.  Variable  varLag  shows  the 
analysis  for  the  output  of  the  variable  lag  filter.  The  results  are  similar  to  that  of  the  example  in 
Section  4. 2. 1.2  and  do  not  require  further  discussion. 
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4.2.3.  Filter  Example  using  a  Lag  Filter 


Figure  19  Lag  Filter  -  Test  Model  Diagram 

Figure  19  shows  ‘lag  filter  quantized’,  a  model  with  a  lag  filter  with  a  fixed  time  constant.  The 
next  subsection  describes  the  transfer  function  and  difference  equation  derivation  for  this  class  of 
filters;  followed  by  the  transfer  function  equations.  Following  subsection  shows  the  analysis 
results  assuming  the  filter  runs  for  72,000  seconds  at  20  Hz  (50  ms  period). 

4.2. 3.1.  Lag  Filter  Transfer  Function 

Figure  20  shows  the  Lag  filter  as  it  appears  as  a  block  in  MATLAB  Simulink.  This  class  of 
filters  has  a  time  constant  tau  that  is  set  for  each  filter  instance.  It  is  not  resettable. 


> 


1 


Tau  s  +1 


lag 


Figure  20  Lag  Filter 

The  following  is  the  transfer  function  derivation  for  this  filter  class: 

Transfer  Function: 


General  form: 

Continuous  domain:  H(s)  —  — —  Discrete  domain: 

V  J  TS+l 

Coefficient  definitions  (all  coefficients  are  constant) 

2t 

fci  =  1  +  —  k2  -  1 

*  S 


H(z)  = 


1+z-1 

fci+/C2z_1 


Difference  equation: 

yn  =  («n  +  «n-l  -  *2yn-l)/fc  1 


(188) 
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4. 2. 3. 2.  Analysis  results  after  72,000  Seconds  (20  hrs) 


System  OperTime  (secs):  72000  Model  Period  i 

ms):  50 Total  Periods:  1440000 

it  Symbol  Name 

Normal 

Range 

Min 

Normal 

Range 

Max 

Range 

Constrai 

nt  Min 
(if  any) 

Range 

Constrai 

nt  Max 

(if  any) 

Data  Type 

Worst  Case 

Numerical  Error 

Bounds 

h_cabin_rate_out_fpm 

0 

225 

-20000 

20000 

Langlndep_Float32 

eB[-+O.0Oe+OOO] 

#  Filters 

Iag3 

11.999944 

235.00006 

Langlndep_Float32 

eB[-+8.55e-010] 

Figure  21  Lag  Filter  -  72000  Seconds 

In  the  above  figure,  hcabinrateoutfpm  has  numerical  error  bounds  of  zero  because  it  is 
processed  through  a  ‘round’ ing  block  .  The  only  other  input  affecting  the  value  at 
h  cabin  rate  out  fpm  is  a  whole  number  (the  RATEQUANTIZER  constant  of  25),  which  can 
be  represented  without  loss.  The  values  for  lag3  above  show  the  analysis  as  applicable  to  the 
output  of  the  lag3  filter.  The  results  are  similar  to  that  of  the  example  in  Section  4. 2. 1.2  and  do 
not  require  further  discussion. 
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4.2.4.  Filter  Example  using  a  Quadratic  Filter 


5witch3  switch* 


Figure  22  Quadratic  Filter  -  Test  Model  Diagram 

Figure  22  shows  ‘sensor  feedback  selection’,  a  model  with  two  quadratic  filters.  The  next 
subsection  describes  the  transfer  function  and  difference  equation  derivation  for  this  class  of 
filters.  Following  subsection  shows  the  analysis  results  assuming  the  filter  runs  for  72,000 
seconds  at  10  Hz  (100  ms  period). 


4.2.4. 1.  Quadratic  Filter  Transfer  Function 


Figure  23  shows  the  Quadratic  filter  as  it  appears  as  a  block  in  MATLAB  Simulink.  This  class  of 
filters  has  six  constants  specified  for  each  instance:  a,b,c,  p,q,r  that  fonn  the  coefficients  of  two 
quadratic  equations,  a  numerator:  asA2  +  bs+  c,  and  a  denominator:  psA2+qs+r.  There  is  also  a 
reset  (IC)  that  resets  the  filter  to  a  pre-specified  value. 


> 

> 


ln  asA2+bs+c 
psA2+qs+r 


> 


quadraticFilter 


Figure  23  Quadratic  Filter 

The  following  is  the  transfer  function  derivation  for  this  filter  class: 


Reset  Mode  (when  IC  =  1):  y  —  u*  cfr  ( c/r  is  the  DC  gain  of  the  filter) 
Filter  Mode  Transfer  Function  (when  IC  —  0): 
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General  form: 


nc^  +  hc+r  J-I-R7  l-i .f  7  2 

Continuous  domain:  //(s)  =  — ; -  Discrete  domain:  H(z )  = - ; - ? 

ps2+qs+r  D+Ez  r+Fz  2 

Coefficient  definitions  (all  coefficients  are  constant) 

A  =  4a  +  2bTs  +  cT2  B  = -2(4a  -  cT2s)  C  =  4a  -  2bTs  +  cT2s 

D  =  Ap  +  2qTs  +  rT2  E  =  -2{Ap-rT 2)  F  =  4p  -  2qTs  +  rT2 

Difference  equation: 

yn  =  (Aun  +  Bun _!  +  Cun_2  -  Eyn_x  -  Fyn_2)/D  (189) 

4.2.4. 2.  Analysis  results  after  72,000  Seconds  (20  hrs) 


System  OperTime  (secs):  72000  Model  Period  (ms):  100  Total  Periods:  720000 

U  Symbol  Name 

Normal 

Range  Min 

Normal 

Range  Max 

Data  Type 

Worst  Case 

Numerical  Error 

Bounds 

DISCHARGED 

-49.99929177 

32.00013311 

Langlndep_Float32 

eB[-+6.58e-003] 

MAX_C 

-49.99929177 

32.00013311 

Langlndep_Float32 

eB[-+6.58e-003] 

max_unfilt_c 

-50 

32 

Langlndep_Fl°at32 

eB[-+3.81e-006] 

FAIL_RECONFIG 

0 

1 

Langlndep_B°°l 

eB[-+O.00e+C00] 

#  Filters 

quadraticFilter 

-25.00008124 

32.00013311 

Langlndep_Float32 

eB[-+3.12e-003] 

quadraticFilterl 

-49.99929177 

-10.00002833 

Langlndep_Float32 

eB[-+6.58e-003] 

Figure  24  Quadratic  Filter  -  72000  Seconds 
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The  range  bounds  of  the  filter  output  variables  in  the  above  figure  are  compliant  with  the 
theoretically  expected  results,  similar  to  that  of  the  example  in  Section  4.2. 1.2,  and  do  not  require 
further  discussion. 

In  the  above  figure,  FAIL  CONFIG  has  numerical  error  bounds  of  zero  because  it  is  a  Boolean. 
DISCHARGEC  and  MAXC  have  the  same  numerical  error  bounds  because  they  come  from 
the  results  of  the  two  quadratic  filters.  MAX  UNFILT  C  has  the  error  bound  of  a  32-bit  floating 
point  representation  because  it  uses  values  straight  from  input  without  any  processing.  As  you 
can  see  the  error  bounds  on  MAX  C  and  DISCHARGE  C  are  the  max  of  the  error  bounds  on 
the  two  quadratic  filters,  since  both  of  them  can  end  up  with  a  value  from  either  filter.  It  is 
remarkable  that  our  static  analysis  techniques  were  able  to  prove  numerical  error  bounds  on  the 
2nd  order  (quadratic)  filter.  The  error  bound  of  [-6.58e-003,  6.58e-003]  on  quadraticFilterl  output 
is  within  the  practical  tolerance  guideline  of  1%  of  the  range.  We  conducted  the  analysis  run  for 
72  seconds,  720  seconds,  7200  seconds,  and  72000  seconds  of  operation.  In  all  cases  the  error 
bounds  is  the  same  ([-6.58e-003,  6.58e-003]);  implying  it  converges  and  stabilizes  over  time; 
thus  supporting  the  theoretical  analysis  expectations. 
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4.2.5.  Filter  Example  using  a  Variable  Washout  Filter 


Figure  25  Variable  Washout  Filter  -  Test  Model  Diagram 

Figure  25  shows  ‘VarWashoutFilter’,  a  model  with  a  variable  washout  filter.  The  next  section 
(subsection  4.2.5. 1)  shows  the  transfer  function  diagram  Figure  26  followed  by  the  transfer 
function  equations.  Subsection  4. 2. 5. 2  shows  the  analysis  results  assuming  the  filter  runs  for 
72,000  seconds  at  10  Hz  (100  ms  period).  Subsection  4.2. 5. 3  shows  the  actual  ‘C’  code  for  the 
filter. 

4.2.5. 1.  Variable  Washout  Filter  Transfer  Function 

Figure  26  shows  the  Variable  Washout  filter  as  it  appears  as  a  block  in  MATLAB  Simulink.  This 
class  of  filters  has  a  variable  tau,  x,  that  determines  the  lag  according  to  the  transfer  function. 
There  are  additional  reset  (IC)  and  reset  value  (ICV)  inputs  that  allow  the  filter  to  be  reset. 

>  In 

>ICV  15  > 

>IC  T5  +  I  ' 

>  7 

warWashout 


ICV 


15 


lc  T5 

7 


Figure  26  Variable  Washout  Filter 

The  following  is  the  transfer  function  derivation  for  this  filter  class: 

Reset  Mode  (when  IC  =  1):  y  —  ICV 
Filter  Mode  Transfer  Function  (when  IC  —  0): 

General  form: 

Continuous  domain:  H(s )  =  Discrete  domain:  H(z )  =  k3+k*z 

v  ts+1  v  '  k1+k2z~1 

Coefficient  definitions  (r  is  given  as  an  input) 
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k-i  —  1  +  — - 

T< 


Difference  equation: 

yn  =  (k3un  +  fc4Mn_!  -  k2yn_x)lk\  (190) 

4.2. 5. 2.  Analysis  results  after  72,000  Seconds  (20  hrs) 


System  OperTime  (secs):  72000  Model  Period  (ms):  100  Total  Periods:  720000 

#  Symbol  Name 

Normal 

Range  Min 

Normal 

Range  Max 

Data  Type 

Worst  Case 

Numerical  Error 

Bounds 

output 

2811.069815 

62188.08073 

Langlndep_Float32 

eB[-+l.lle+001] 

#  Filters 

varWashout 

937.0232717 

20729.36024 

Langlndep_Float32 

eB[-+3.71e+000] 

Figure  27  Variable  Washout  Filter  -72000  Seconds 

In  the  above  figure,  you  can  see  that  output  has  an  error  bound  of  3  times  the  error  bound  of 
varWashout,  because  they  are  separated  by  a  gain  of  3. 
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4.3.  Summary  of  Accumulator  Test  Suite  Results 


The  development  of  the  second  year  analysis  of  Integrating  Accumulators  was  guided  by,  tuned, 
and  ultimately  tested  for  round  trip  integration  by  a  suite  of  test  cases  designed  by  the  Honeywell 
team.  These  consisted  of  Simulink  models  of  representative  accumulator  applications,  their 
generated  C  implementations  and  their  generated  JSON  files  describing  model-level  properties 
of  the  generated  code  to  be  exploited  by  CodeHawk.  This  section  summarizes  the  results  of  these 
test  cases. 

Motivation  and  Significance  of  the  Results 

The  Patriot  Missile  system  suffered  an  infamous  failure  [5]  when  its  computation  of  time  drifted 
over  the  course  of  two  days  of  continuous  operation.  The  computations  using  the  24-bit  floating 
point  unit  introduced  error  in  the  calculation  of  elapsed  time  accumulated  since  system  start. 
After  8  hours  of  operations,  this  resulted  in  a  20  percent  error  in  the  "range  gate"  -  the  expected 
position  of  the  target.  After  20  hours,  the  error  was  as  much  as  50  percent,  resulting  the  Patriot 
missile’s  failure  to  intercept  an  incoming  Scud  missile.  Figure  28  shows  the  increase  in  absolute 
error  with  elapsed  time. 


Hours 

Seconds 

Calculated  Time  (sec) 

Inaccuracy  (sec) 

Approx,  shift  in 
Range  Gate 
(meters) 

0 

0 

0 

0 

0 

1 

3600 

3599.9966 

.  0034 

7 

8 

28800 

8799.9725 

.  0251 

55 

20  (a) 

72000 

71999.9313 

.0687 

137 

48 

172800 

172799.8352 

.1648 

330 

72 

259200 

259199.7528 

.2472 

494 

100  (b) 

360000 

359999.6667 

.3433 

687 

Figure  28.  The  error  in  accumulated  time  due  to  floating  point  multiplication  and  the  resulting  distance  by  which  the 
computed  range  gate  was  off 


Using  the  combination  of  model-level  and  code-level  analysis,  we  have  addressed  the  problem  of 
automated  detection  and  analysis  of  such  constructions.  HiLiTE  performance  an  analysis  of 
cyclic  path  (feedback  loop)  invariants  in  the  Simulink  model  to  detect  a  counter  type  of 
invariants  and  provide  this  information  to  CodeHawk  in  terms  of  the  accumulator  variable  and 
nature  of  the  invariant.  CodeHawk  then  applies  novel  code  analysis  techniques  to  derive  precise 
closed-form  mathematical  bounds  for  the  both  the  range  and  numerical  error  as  a  function  of 
time.  The  results  presented  in  the  following  subsections  clearly  show  how  the  numerical  error 
grows  with  system  operational  time. 
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4.3.1.  Accumulator  Example  using  a  Fixed  Integer  Increment 


simpleDelay  c=double(0.l) 


{IC=0) 


time  (a)  [11 


inport  (f)  [1] 
10  to  20 


inportl  (f)  [2] 
0  to  5 


► 

► 

suml 


^  outport  (a)  [2] 


Figure  29.  Fixed  Increment  Accumulator  -  abstract  diagram  for  the  Patriot  Missile  Bug 

Figure  29  shows  a  model  with  an  accumulator  with  a  fixed  integer  increment  which  abstracts  the 
essential  elements  of  the  Patriot  Missile  bug.  At  each  periodic  execution  frame  (at  10  Hz  rate) 
the  counter  (output  of  sum  block)  counts  by  1 .  The  numerical  error  occurs  when  the  counter 
value  is  multiplied  by  0. 1  to  get  the  value  of  elapsed  time;  note  that  0. 1  is  not  accurately 
representable  in  any  binary-based  floating  point  notation.  This  numerical  error  is  shown  in  the 
results. 

The  next  subsection  shows  the  analysis  results  assuming  the  filter  runs  for  72,  720,  7200  and 
72,000  seconds  at  10  Hz  (100  ms  period).  The  JSON  files  exchanged  between  HiLiTE  and 
CodeHawk  and  the  actual  ‘C’  code. 

4.3. 1.1.  Analysis  results  after  different  periods  (72,  720,  7200,  72000  seconds) 


System  OperTime  (secs):  72  Model  Period  (ms) 

:  lOO  Total  Periods:  720 

Accumulator  pattern  detected  - 

Path  Invariant  Analysis: 

Block  :  sum;  Variable  :  y(k)_l;  Invariant: 

Guard:  Assignment:  y(k+l)_l=l+y(k)_l; 

#  Symbol  Name 

Normal  Range 
Min 

Normal  Range 

Max 

Data  Type 

Worst  Case 

Numerical  Error 

Bounds 

time 

O.l 

72.00000107 

Langlndep_Float32 

eB[-+4.30e-005] 

outport 

lO 

25 

Langlndep_Float32 

eB[-+2.38e-006] 

Figure  30  Fixed  Increment  Accumulator  -  72  Seconds 
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System  OperTime  (secs):  720  Model  Period  (ms):  lOO  Total  Periods: 

7200 

Accumulator  pattern  detected  - 

Path  Invariant  Analysis: 

Block  :  sum;  Variable  :  y(k)_l;  Invariant: 

Guard:  Assignment:  y(k+l)_l=l+y(k)_l; 

#  Symbol  Name 

Normal  Range 
Min 

Normal  Range 

Max 

Data  Type 

Worst  Case 

Numerical  Error 

Bounds 

time 

O.l 

720.0000107 

Langlndep_Float32 

eB[-+4.29e-004] 

outport 

lO 

25 

Langlndep_Float32 

eB[-+2.38e-006] 

Figure  31  Fixed  Increment  Accumulator  -  720  Seconds 


System  OperTime  (secs):  7200  Model  Period  (ms):  lOO  Total  Periods:  72000 

Accumulator  pattern  detected  -  Path  Invariant  Analysis: 

Block  :  sum;  Variable  :  y(k)_l;  Invariant: 

Guard:  Assignment:  y(k+l)_l=l+y(k)_l; 

U  Symbol  Name 

Normal  Range 
Min 

Normal  Range 

Max 

Data  Type 

Worst  Case 

Numerical  Error 

Bounds 

time 

0.1 

7200.000107 

Langlndep_Float32 

eB[--H4.31e-003] 

outport 

lO 

25 

Langlndep_Float32 

eB[-+2.38e-006] 

Figure  32  Fixed  Increment  Accumulator  -  7200  Seconds 


System  OperTime  (secs):  72000  Model  Period  (ms):  lOO  Total  Periods:  720000 

Accumulator  pattern  detected  - 

Path  Invariant  Analysis: 

Block  :  sum;  Variable  :  y(k)_l;  Invariant: 

Guard:  Assignment:  y( k+l)_l=l+y(k)_l; 

#  Symbol  Name 

Normal  Range 
Min 

Normal  Range 

Max 

Data  Type 

Worst  Case 

Numerical  Error 

Bounds 

time 

0.1 

72000.00107 

Langlndep_Float32 

eB[--H4.48e-002] 

outport 

lO 

25 

Langlndep_Float32 

eB[-+2.38e-006] 

Figure  33  Fixed  Increment  Accumulator  -  72000  Seconds 


In  the  above  figure,  the  range  and  error  bounds  for  the  time  variable  are  provide  for  different 
system  operational  times;  both  grow  as  a  linear  function  of  time  as  theoretically  expected,  (note 
that  variable  outport  in  this  example  is  unrelated  to  the  accumulator  and  thus  its  error  bound 
remains  fixed  with  the  operational  time  ). 
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4.3.2.  Accumulator  Example  using  Variable  increments 


increment  (f)  [1] 


max  allow  range:  [0.1, 0.2] 


Figure  34  Variable  Increment  Accumulator  -  Test  Model  Diagram 

The  figure  above  shows  ‘CounterWithoutLimitsWithExtemallnputs’,  a  model  with  a  floating¬ 
point  accumulator  with  a  variable  increment  with  the  input  range  constrained  to  [0.1,  0.2]  .  The 
next  subsection  shows  the  analysis  results  assuming  the  filter  runs  for  72,  720,  7200  72,000 
seconds  at  20  Hz  (50  ms  period).  This  example  was  created  to  ensure  that  both  HiLiTE  can 
CodeHawk  can  properly  identify  and  analyze  accumulators  with  floating-point  variable 
increment.  The  results  are  consistent  with  the  theoretical  analysis. 


4.3. 2.1.  Analysis  results  after  different  periods  (72,  720,  7200,  72000  seconds) 

System  OperTime  (secs):  72  Model  Period  (ms):  SOTotal  Periods:  1440 


Accumulator  pattern  detected  - 

Path  Invariant  Analysis: 

Block  :  sum;  Variable  :  y(k)_l;  Invariant: 

Guard:  Assignment:  y(k+l)_l=u(k+l)_0+y(k)_l; 

ff  Symbol  Name 

Normal  Range 
Min 

Normal  Range 

Max 

Data  Type 

Worst  Case 

Numerical  Error 

Bounds 

output 

0.32 

921.5999794 

Langlndep_Float32 

eB[-+6.60e-004] 

#  Accumulators 

sum 

O.l 

288.0000043 

Langlndep_Float32 

eB[-+1.37e-004] 

Figure  35  Variable  Increment  Accumulator  -  72  Seconds 
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System  OperTime  (secs):  720  Model  Period  (ms):  50  Total  Periods:  14400 

Accumulator  pattern  detected  - 

Path  Invariant  Analysis: 

Block  :  sum;  Variable  :  y(k)_l;  Invariant: 

Guard:  Assignment:  y(k+l)_l=u(k+l)_0+y(k)_l; 

#  Symbol  Name 

Normal  Range 

Min 

Normal  Range 

Max 

Data  Type 

Worst  Case 

Numerical  Error 

Bounds 

output 

0.32 

9215.999794 

Langlndep_Float32 

eB[-+6.60e-003] 

ft  Accumulators 

sum 

0.1 

2880.000043 

Langlndep_Float32 

eB[-+1.37e-003] 

Figure  36  Variable  Increment  Accumulator  -  720  Seconds 


System  OperTime  (secs):  7200  Model  Period  (ms):50Total  Periods:  144000 

Accumulator  pattern  detected  - 

Path  Invariant  Analysis: 

Block  :  sum;  Variable  :  y(k)_l;  Invariant: 

Guard:  Assignment:  v(k+l)  l=u(k+l)  0+y(k)  1; 

#  Symbol  Name 

Normal  Range 

Min 

Normal  Range 

Max 

Data  Type 

Worst  Case 

Numerical  Error 

Bounds 

output 

0.32 

92159.99794 

Langlndep  Float32 

eB[-+6.65e-002] 

#  Accumulators 

sum 

0.1 

28800.00043 

Langlndep_Float32 

eB[-+1.39e-002] 

Figure  37  Variable  Increment  Accumulator  -  7200  Seconds 
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System  OperTime  (secs):  72000  Model  Period  (ms):  50  Total  Periods: 

1440000 

Accumulator  pattern  detected  -  Path  Invariant  Analysis: 

Block  :  sum;  Variable  :  y(k)_l;  Invariant: 

Guard:  Assignment:  y(k+l)_l=u(k-KL)_0+y(k)_l; 

#  Symbol  Name 

Normal  Range 

Min 

Normal  Range 

Max 

Data  Type 

Worst  Case 

Numerical  Error 

Bounds 

output 

0.32 

921599.9794 

Langlndep_Float32 

eB[-+7.19e-001] 

#  Accumulators 

sum 

0.1 

288000.0043 

Langlndep_Float32 

eB[-+1.50e-001] 

Figure  38  Variable  Increment  Accumulator  -  72000  Seconds 


5.  CONCLUSIONS 

We  have  demonstrated  a  combined  model-level  and  code-level  analysis  approach  to  prove 
precise,  yet  conservative  bounds  for  variable  ranges  and  numerical  error  on  two  classes  of 
numerical  algorithms:  linear  digital  filters  and  accumulators.  Our  novel  analysis  approach  uses 
the  semantics  of  control  algorithms  available  in  the  Simulink  models  to  identify  and  classify  the 
properties  of  filter  blocks  and  accumulator  patterns.  The  HiLiTE  tool  introduced  a  novel 
technique  to  analyze  cyclic  path  (feedback  loop)  expressions  in  models  to  automatically  detect 
patterns  and  reduce  them  to  temporal  invariants  in  terms  of  a  single  accumulator  variable  that 
could  then  be  provided  to  CodeHawk.  We  defined  extensible  formats  for  representing  variable 
ranges,  properties  of  filters  and  accumulators,  functions  for  error  bounds.  These  were  used  in 
developing  a  flexible  JSON  interchange  infrastructure  in  the  tool  chain  that  allowed  the 
automated  round-trip  analysis  and  presentation  of  the  results  to  the  user. 

At  the  code-level  analysis,  Anca  Browne  discovered  a  technique,  based  on  classical  discrete 
mathematics,  for  computing  the  exact  (real)  bounds  on  linear  discrete  filters,  and  applied  it  to 
first-order  and  second-order  filters  supplied  by  Honeywell.  The  time  to  analyze  a  filter  is 
negligible.  This  contrasts  with  analyses  performed  for  Airbus  that  take  on  the  order  of  20  hours 
to  get  less  precise  results. 

Extending  this  technique  to  floating  point  error  analysis,  she  discovered  how  to  generate  a 
closed- form  solution  for  the  error  as  a  function  of  the  number  of  iterations  of  the  recurrence  that 
are  performed.  Again,  this  result  was  applied  to  first-order  and  second-order  filters  supplied  by 
Honeywell.  The  result  of  the  error  analysis  were  sufficiently  precise  and  valid  even  for  second- 
order  filters  the  numerical  error  stabilized  within  practically  acceptable  tolerance  bounds  even 
for  very  large  periods  of  run  time.  These  precise  bounds  prove  the  essential  characteristics  of 
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variable  ranges  and  numerical  error  in  the  control  algorithms,  as  discussed  in  detail  in  Section 
4.2. 


These  techniques  were  slightly  generalized  to  analyze  a  different  class  of  numerical  algorithms 
used  in  control  systems,  namely  accumulators.  Accumulators,  such  as  counters,  are  used  in 
many  control  systems  -  a  prime  example  is  the  Patriot  Missile  system  where  the  accumulated 
numerical  error  resulted  in  an  infamous  failure.  Again,  our  techniques  allow  almost 
instantaneous  generation  of  an  analysis  function  that  gives  (1)  accumulator  results  as  a  function 
of  the  number  of  iterations  performed,  and  (2)  accumulated  floating-point  error  as  a  function  of 
the  number  of  iterations  performed. 

As  avionics  software  becomes  increasingly  complex,  the  need  for  automated  analysis  becomes 
critical  to  verify  correct  algorithm  behavior  and  absence  of  potential  failures.  The  current  static 
analysis  tools  and  proposed  techniques  may  handle  discrete  types  of  algorithms  well,  but  fall 
short  of  providing  a  precise  analysis  for  numerically  intensive  algorithms.  This  can  result  in  large 
range  bounds  (or  no  error  bounds)  provided  by  the  tools  which  then  leave  it  to  the  developer  to 
analyze  the  problem  by  manual  reviews  or  empirical  means.  This  is  already  proving  to  be  a 
problem  in  current  avionics  system  when  deploying  large  amounts  of  software  with  frequent 
configuration  changes.  For  example,  developers  may  need  to  do  extensive  analysis  to  ascertain 
that  the  new  set  of  values  of  filter  parameters  is  not  degenerate  and  will  indeed  make  the  range 
and  error  bounds  converge. 

We  believe  this  work  is  a  significant  step  towards  improving  the  automation  and  quality  of 
results  provided  by  static  analysis  techniques.  Our  novel  analysis  techniques  can  automatically 
analyze  the  code  implementations  of  a  variety  of  first  and  second  order  linear  filter  and  derive 
not  only  precise  range  bounds  but  also  numerical  error  bounds.  This  allows  quick  verification  of 
normal  algorithm  behavior  as  well  as  detection  of  erroneous  combinations  filter  parameters  (time 
constant  values)  that  can  make  bounds  diverge,  and  defective  implementations  that  can  give  rise 
to  large  numerical  error.  Our  techniques  for  the  analysis  of  accumulator  patterns  would  have 
detected  the  bug  that  led  to  Patriot  Missile  failure,  and  may  hopefully  prevent  such  faulty 
software  in  the  future. 
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7.  LIST  OF  ACRONYMS 


COM 

Common  Object  Model 

HiLiTE 

Honeywell  Integrated  Lifecycle  Tool  Environment 

JSON 

JavaScript  Object  Notation 

LDF 

Linear  Digital  Filters 

MBD 

Model-Based  Development 

Ocaml 

0  Collaborative  Application  Markup  Language 

PI 

Proportional-Integrator 

SANA 

Static  Analysis  of  Numerical  Algorithms 

SoW 

Statement  of  Work 

WCEB 

Worst  Case  Error  Bound 

Approved  for  Public  Release;  Distribution  Unlimited. 

67 


